Data

Morphological data

Morphological measurements (eye diameter and body length) were taken from a mixture of fresh and preserved specimens. In total, 459 sergestid specimens were sampled from 16 species belonging to 12 genera and two subgroups (Sergestes and Sergia). Measurements were recorded to ±0.01mm, but we rounded all to the nearest ±0.1mm as we were not confident in the accuracy of such fine measurements taken with calipers on a rocking ship.

# Full specimen morphology dataset  ------

#import data
specimens_raw <- data.frame(read.csv("../Data/specimen_data.csv", header=TRUE, na.strings=c("", "NA", " ", stringsAsFactors = FALSE)))

#tidy data
specimens <- specimens_raw %>%
  mutate(Genus = gsub('\\s+', '', Genus)) %>%
  mutate(Family = str_to_title(Family)) %>%
  mutate(Subgroup = str_to_title(Subgroup)) %>%
  mutate(Genus = str_to_title(Genus)) %>%
  mutate(genus_species = as.factor(paste(Genus, Species, sep = "_"))) %>%
  mutate(pres_bin = recode_factor(Preservation, "ethanol" = "preserved", "paraformaldehyde" = "preserved", "fresh" = "fresh")) %>%
  mutate_if(is.numeric,round, 1)


# Species traits  ------

#import data
traits_raw <- data.frame(read.csv("../Data/species_data.csv", header=TRUE, na.strings=c("", "NA", " ", stringsAsFactors = FALSE)))

#tidy data
traits <- traits_raw %>%
  mutate(Genus = gsub('\\s+', '', Genus)) %>%
  mutate(genus_species = as.factor(paste(Genus, Species, sep = "_"))) %>%
  mutate(Range = Day_Avg - Night_Avg)

# Show sampling -----

#compile number of shrimp sampled by species
counts <-ddply(specimens, .(specimens$Subgroup, specimens$Genus, specimens$Species), nrow) 

#create table column names
names(counts) <- c("Subgroup", "Genus","Species","Sampled") 

#generate scrolling table in RMarkdown
kable(counts[ , c("Subgroup", "Genus","Species","Sampled")], caption = "Species and sampling effort for morphological data from family Sergestidae.") %>% 
                      kable_styling(full_width = F) %>% 
                      collapse_rows(columns = 1, valign = "top") %>%
                      scroll_box(height = "400px") 
Species and sampling effort for morphological data from family Sergestidae.
Subgroup Genus Species Sampled
Sergestes Allosergestes pectinatus 14
Allosergestes sargassi 38
Deosergestes corniculum 18
Deosergestes henseni 58
Eusergestes arcticus 5
Neosergestes edwardsii 7
Parasergestes armatus 30
Parasergestes vigilax 17
Sergestes atlanticus 15
Sergia Challengerosergia hansjacobi 43
Challengerosergia talismani 32
Gardinerosergia splendens 42
Phorcosergia grandis 51
Robustasergia robusta 23
Robustosergia regalis 38
Sergia tenuiremis 28

Molecular phylogeny

Below is the sergestid molecular phylogeny generated by Charles Golightly. We used ape v.5.4-1 to prune the tree to the 16 species in our dataset and make the tree ultrametric (using the chronopl function with lambda set to 1).

#Import phylogeny
tree_orig <- read.tree(file = "../Data/sertestid_tree_rooted")

#Make copy of tree to manipulate
tree_edit <- tree_orig

#Remove locality tags from phylo tip labels 
n <- 3 #which _ to end pattern at
pat <- paste0('^([^_]+(?:_[^_]+){',n-1,'}).*') #pattern
tree_edit$tip.label <- sub(pat, '\\1', tree_edit$tip.label) #remove tags from tips

#Make list of taxa to drop (in tree but not in dataset)
drops <- setdiff(tree_edit$tip.label, as.character(traits$Binomial))

#Drop unwanted tips from phylogeny
tree.pruned <- drop.tip(phy = tree_edit, tip = drops) 

#Check and reorder tree
is.rooted(tree.pruned) 
is.binary.tree(tree.pruned) 
tree <- ladderize(tree.pruned) 

#Make row names of the species the phylogeny tip labels
rownames(traits) <- traits$Binomial

#Reorder species data to match phylogeny order
traits <-ReorderData(tree, traits, taxa.names="row names")

#check that phylogeny tips and data match exactly (if they match will return "OK")
name.check(tree, traits)

#make tree ultrametric
tree <- chronopl(tree, 1)
#plot tree of sampled species with original tip labels
plot.phylo(tree, #phylogeny to plot
           type = "phylogram", #how to visualize phylogeny
           show.tip.label = TRUE, #whether to show tip labels/species names
           cex = 1.0, #text size
           no.margin = TRUE, 
           use.edge.length = TRUE,
           edge.width = 2, #thickness of phylogeny branches
           label.offset = 0.1) #how far away from tips to put labels

Static eye-body allometry in Sergestidae

Eye-body allometry in each species

We first examine static post-metamorphosis eye diameter-body length allometry across 16 species of deep-sea sergestids. Allometric relationships were fit using standardised major axis regression via the smatr v.3.4.8 package. For each species, standard model checks, model fits, and an interactive plot of data and fits are shown. Hover over the plots to see individual data values. Points are colored by the source of the sample (fresh collection or preserved).

Note that for 3 of the species, a large/unreasonable outlier was present. We identified these by eye and they have been removed from the data subsets and analyses, but the outlier value is noted under each species header for reference, and a second plot including the outlier as a red diamond is included for transparency (these can go into the supplemental information if we want).

#set colors for fresh and preserved samples
col_pres <- c("ethanol" = "#a6cee3",
              "fresh" = "#b2df8a",
              "paraformaldehyde" = "#1f78b4")

Allosergestes pectinatus

# Subset data -----
apec <- specimens %>% filter(genus_species == "Allosergestes_pectinatus")

# Fit SMA model ------
sma_apec <- sma(formula = Eye_Diameter ~ Body_Length, 
            data = apec, 
            log = "xy", #sets both x and y variables as logged
            method="SMA", #defines SMA as model
            alpha = 0.05)

#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_apec, which = "default",type = "o") 
plot(sma_apec, which = "residual",type = "o")
plot(sma_apec, which = "qq", type = "o")
par(mfrow = c(1,1)) 

#view fits
summary(sma_apec)
## Call: sma(formula = Eye_Diameter ~ Body_Length, data = apec, log = "xy", 
##     method = "SMA", alpha = 0.05) 
## 
## Fit using Standardized Major Axis 
## 
## These variables were log-transformed before fitting: xy 
## 
## Confidence intervals (CI) are at 95%
## 
## ------------------------------------------------------------
## Coefficients:
##             elevation    slope
## estimate    -2.812398 1.940053
## lower limit -3.491769 1.485663
## upper limit -2.133027 2.533418
## 
## H0 : variables uncorrelated
## R-squared : 0.8156794 
## P-value : 9.643e-06
#save coefficients of fits as object
cc_sma_apec <- data.frame(coef(sma_apec))

# Make plot ------
plot_apec <- ggplot(apec, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) + 
  geom_point(size = 2, alpha = 0.8) + 
  scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
  scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
  scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
  ggtitle("Allosergestes pectinatus") + 
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
  geom_abline(data = cc_sma_apec, aes(intercept = cc_sma_apec[1,1], slope = cc_sma_apec[2,1])) +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_apec[2,1], digits = 2), sep = " "), size = 3, col = "black") +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_apec[1,1], digits = 2), sep = " "), size = 3, col = "black")

ggplotly(plot_apec)

Allosergestes sargassi

Note: this species has a big outlier that was removed for analyses and plots (eye 0.7, body 13.0).

# Subset data -----
#note: outlier removed
asar <- specimens %>% 
  filter(genus_species == "Allosergestes_sargassi") %>%
  filter(!(Eye_Diameter == round(0.67, 1) & Body_Length == round(12.97,1)))

# Fit SMA model ------
sma_asar <- sma(formula = Eye_Diameter ~ Body_Length, 
            data = asar, 
            log = "xy", #sets both x and y variables as logged
            method="SMA", #defines SMA as model
            alpha = 0.05)

#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_asar, which = "default",type = "o") 
plot(sma_asar, which = "residual",type = "o")
plot(sma_asar, which = "qq", type = "o")
par(mfrow = c(1,1)) 

#view fits
summary(sma_asar)
## Call: sma(formula = Eye_Diameter ~ Body_Length, data = asar, log = "xy", 
##     method = "SMA", alpha = 0.05) 
## 
## Fit using Standardized Major Axis 
## 
## These variables were log-transformed before fitting: xy 
## 
## Confidence intervals (CI) are at 95%
## 
## ------------------------------------------------------------
## Coefficients:
##             elevation    slope
## estimate    -2.728530 1.844597
## lower limit -3.402995 1.419494
## upper limit -2.054065 2.397008
## 
## H0 : variables uncorrelated
## R-squared : 0.4037707 
## P-value : 2.3852e-05
#save coefficients of fits as object
cc_sma_asar <- data.frame(coef(sma_asar))

# Make plot ------
plot_asar <- ggplot(asar, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) + 
  geom_point(size = 2, alpha = 0.8) + 
  scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
  scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
  scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
  ggtitle("Allosergestes sargassi") + 
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
  geom_abline(data = cc_sma_asar, aes(intercept = cc_sma_asar[1,1], slope = cc_sma_asar[2,1])) +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_asar[2,1], digits = 2), sep = " "), size = 3, col = "black") +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_asar[1,1], digits = 2), sep = " "), size = 3, col = "black")

#plot with outlier shown
plot_asar_out <- plot_asar +
  geom_point(aes(y=round(0.67, 1), x=round(12.97,1)), colour="red", pch = 18) +
  ggtitle("Allosergestes sargassi + outlier")
  
#plot
ggplotly(plot_asar)
#plot with outlier
ggplotly(plot_asar_out)

Challengerosergia hansjacobi

# Subset data -----
chan <- specimens %>% filter(genus_species == "Challengerosergia_hansjacobi")

# Fit SMA model ------
sma_chan <- sma(formula = Eye_Diameter ~ Body_Length, 
            data = chan, 
            log = "xy", #sets both x and y variables as logged
            method="SMA", #defines SMA as model
            alpha = 0.05)

#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_chan, which = "default",type = "o") 
plot(sma_chan, which = "residual",type = "o")
plot(sma_chan, which = "qq", type = "o")
par(mfrow = c(1,1)) 

#view fits
summary(sma_chan)
## Call: sma(formula = Eye_Diameter ~ Body_Length, data = chan, log = "xy", 
##     method = "SMA", alpha = 0.05) 
## 
## Fit using Standardized Major Axis 
## 
## These variables were log-transformed before fitting: xy 
## 
## Confidence intervals (CI) are at 95%
## 
## ------------------------------------------------------------
## Coefficients:
##             elevation    slope
## estimate    -1.528299 1.047036
## lower limit -1.913161 0.824145
## upper limit -1.143436 1.330207
## 
## H0 : variables uncorrelated
## R-squared : 0.4129122 
## P-value : 3.3885e-06
#save coefficients of fits as object
cc_sma_chan <- data.frame(coef(sma_chan))

# Make plot ------
plot_chan <- ggplot(chan, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) + 
  geom_point(size = 2, alpha = 0.8) + 
  scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
  scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
  scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
  ggtitle("Challengerosergia hansjacobi") + 
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
  geom_abline(data = cc_sma_chan, aes(intercept = cc_sma_chan[1,1], slope = cc_sma_chan[2,1])) +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_chan[2,1], digits = 2), sep = " "), size = 3, col = "black") +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_chan[1,1], digits = 2), sep = " "), size = 3, col = "black")

ggplotly(plot_chan)

Challengerosergia talismani

# Subset data -----
ctal <- specimens %>% filter(genus_species == "Challengerosergia_talismani")

# Fit SMA model ------
sma_ctal <- sma(formula = Eye_Diameter ~ Body_Length, 
            data = ctal, 
            log = "xy", #sets both x and y variables as logged
            method="SMA", #defines SMA as model
            alpha = 0.05)

#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_ctal, which = "default",type = "o") 
plot(sma_ctal, which = "residual",type = "o")
plot(sma_ctal, which = "qq", type = "o")
par(mfrow = c(1,1)) 

#view fits
summary(sma_ctal)
## Call: sma(formula = Eye_Diameter ~ Body_Length, data = ctal, log = "xy", 
##     method = "SMA", alpha = 0.05) 
## 
## Fit using Standardized Major Axis 
## 
## These variables were log-transformed before fitting: xy 
## 
## Confidence intervals (CI) are at 95%
## 
## ------------------------------------------------------------
## Coefficients:
##             elevation     slope
## estimate    -1.578063 1.0389917
## lower limit -1.840235 0.8845412
## upper limit -1.315891 1.2204109
## 
## H0 : variables uncorrelated
## R-squared : 0.8120894 
## P-value : 2.0471e-12
#save coefficients of fits as object
cc_sma_ctal <- data.frame(coef(sma_ctal))

# Make plot ------
plot_ctal <- ggplot(ctal, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) + 
  geom_point(size = 2, alpha = 0.8) + 
  scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
  scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
  scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
  ggtitle("Challengerosergia talismani") + 
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
  geom_abline(data = cc_sma_ctal, aes(intercept = cc_sma_ctal[1,1], slope = cc_sma_ctal[2,1])) +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_ctal[2,1], digits = 2), sep = " "), size = 3, col = "black") +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_ctal[1,1], digits = 2), sep = " "), size = 3, col = "black")

ggplotly(plot_ctal)

Deosergestes corniculum

# Subset data -----
dcor <- specimens %>% filter(genus_species == "Deosergestes_corniculum")

# Fit SMA model ------
sma_dcor <- sma(formula = Eye_Diameter ~ Body_Length, 
            data = dcor, 
            log = "xy", #sets both x and y variables as logged
            method="SMA", #defines SMA as model
            alpha = 0.05)

#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_dcor, which = "default",type = "o") 
plot(sma_dcor, which = "residual",type = "o")
plot(sma_dcor, which = "qq", type = "o")
par(mfrow = c(1,1)) 

#view fits
summary(sma_dcor)
## Call: sma(formula = Eye_Diameter ~ Body_Length, data = dcor, log = "xy", 
##     method = "SMA", alpha = 0.05) 
## 
## Fit using Standardized Major Axis 
## 
## These variables were log-transformed before fitting: xy 
## 
## Confidence intervals (CI) are at 95%
## 
## ------------------------------------------------------------
## Coefficients:
##             elevation     slope
## estimate    -1.432046 0.8653057
## lower limit -1.737583 0.7020666
## upper limit -1.126509 1.0664999
## 
## H0 : variables uncorrelated
## R-squared : 0.8421211 
## P-value : 8.1768e-08
#save coefficients of fits as object
cc_sma_dcor <- data.frame(coef(sma_dcor))

# Make plot ------
plot_dcor <- ggplot(dcor, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) + 
  geom_point(size = 2, alpha = 0.8) + 
  scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
  scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
  scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
  ggtitle("Deosergestes corniculum") + 
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
  geom_abline(data = cc_sma_dcor, aes(intercept = cc_sma_dcor[1,1], slope = cc_sma_dcor[2,1])) +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_dcor[2,1], digits = 2), sep = " "), size = 3, col = "black") +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_dcor[1,1], digits = 2), sep = " "), size = 3, col = "black")

ggplotly(plot_dcor)

Deosergestes henseni

# Subset data -----
dhen <- specimens %>% filter(genus_species == "Deosergestes_henseni")

# Fit SMA model ------
sma_dhen <- sma(formula = Eye_Diameter ~ Body_Length, 
            data = dhen, 
            log = "xy", #sets both x and y variables as logged
            method="SMA", #defines SMA as model
            alpha = 0.05)

#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_dhen, which = "default",type = "o") 
plot(sma_dhen, which = "residual",type = "o")
plot(sma_dhen, which = "qq", type = "o")
par(mfrow = c(1,1)) 

#view fits
summary(sma_dhen)
## Call: sma(formula = Eye_Diameter ~ Body_Length, data = dhen, log = "xy", 
##     method = "SMA", alpha = 0.05) 
## 
## Fit using Standardized Major Axis 
## 
## These variables were log-transformed before fitting: xy 
## 
## Confidence intervals (CI) are at 95%
## 
## ------------------------------------------------------------
## Coefficients:
##             elevation     slope
## estimate    -1.306994 0.8238705
## lower limit -1.593941 0.6637531
## upper limit -1.020047 1.0226131
## 
## H0 : variables uncorrelated
## R-squared : 0.3380975 
## P-value : 1.6982e-06
#save coefficients of fits as object
cc_sma_dhen <- data.frame(coef(sma_dhen))

# Make plot ------
plot_dhen <- ggplot(dhen, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) + 
  geom_point(size = 2, alpha = 0.8) + 
  scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
  scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
  scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
  ggtitle("Deosergestes henseni") + 
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
  geom_abline(data = cc_sma_dhen, aes(intercept = cc_sma_dhen[1,1], slope = cc_sma_dhen[2,1])) +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_dhen[2,1], digits = 2), sep = " "), size = 3, col = "black") +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_dhen[1,1], digits = 2), sep = " "), size = 3, col = "black")

ggplotly(plot_dhen)

Eusergestes arcticus

# Subset data -----
euar <- specimens %>% filter(genus_species == "Eusergestes_arcticus")

# Fit SMA model ------
sma_euar <- sma(formula = Eye_Diameter ~ Body_Length, 
            data = euar, 
            log = "xy", #sets both x and y variables as logged
            method="SMA", #defines SMA as model
            alpha = 0.05)

#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_euar, which = "default",type = "o") 
plot(sma_euar, which = "residual",type = "o")
plot(sma_euar, which = "qq", type = "o")
par(mfrow = c(1,1)) 

#view fits
summary(sma_euar)
## Call: sma(formula = Eye_Diameter ~ Body_Length, data = euar, log = "xy", 
##     method = "SMA", alpha = 0.05) 
## 
## Fit using Standardized Major Axis 
## 
## These variables were log-transformed before fitting: xy 
## 
## Confidence intervals (CI) are at 95%
## 
## ------------------------------------------------------------
## Coefficients:
##              elevation     slope
## estimate    -2.1539322 1.3306114
## lower limit -3.8635100 0.6694485
## upper limit -0.4443543 2.6447543
## 
## H0 : variables uncorrelated
## R-squared : 0.8368056 
## P-value : 0.029484
#save coefficients of fits as object
cc_sma_euar <- data.frame(coef(sma_euar))

# Make plot ------
plot_euar <- ggplot(euar, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) + 
  geom_point(size = 2, alpha = 0.8) + 
  scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
  scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
  scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
  ggtitle("Eusergestes arcticus") + 
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
  geom_abline(data = cc_sma_euar, aes(intercept = cc_sma_euar[1,1], slope = cc_sma_euar[2,1])) +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_euar[2,1], digits = 2), sep = " "), size = 3, col = "black") +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_euar[1,1], digits = 2), sep = " "), size = 3, col = "black")

ggplotly(plot_euar)

Gardinerosergia splendens

Note: this species has a big outlier that was removed for analyses and plots (eye 1.7, body 22.0).

# Subset data -----
gspl <- specimens %>% 
  filter(genus_species == "Gardinerosergia_splendens") %>%
  filter(!(Eye_Diameter == round(1.71, 1) & Body_Length == round(21.98,1)))

# Fit SMA model ------
sma_gspl <- sma(formula = Eye_Diameter ~ Body_Length, 
            data = gspl, 
            log = "xy", #sets both x and y variables as logged
            method="SMA", #defines SMA as model
            alpha = 0.05)

#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_gspl, which = "default",type = "o") 
plot(sma_gspl, which = "residual",type = "o")
plot(sma_gspl, which = "qq", type = "o")
par(mfrow = c(1,1)) 

#view fits
summary(sma_gspl)
## Call: sma(formula = Eye_Diameter ~ Body_Length, data = gspl, log = "xy", 
##     method = "SMA", alpha = 0.05) 
## 
## Fit using Standardized Major Axis 
## 
## These variables were log-transformed before fitting: xy 
## 
## Confidence intervals (CI) are at 95%
## 
## ------------------------------------------------------------
## Coefficients:
##              elevation     slope
## estimate    -1.2253950 0.8523703
## lower limit -1.5326713 0.6732695
## upper limit -0.9181186 1.0791149
## 
## H0 : variables uncorrelated
## R-squared : 0.4597308 
## P-value : 1.1132e-06
#save coefficients of fits as object
cc_sma_gspl <- data.frame(coef(sma_gspl))

# Make plot ------
plot_gspl <- ggplot(gspl, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) + 
  geom_point(size = 2, alpha = 0.8) + 
  scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
  scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
  scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
  ggtitle("Gardinerosergia splendens") + 
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
  geom_abline(data = cc_sma_gspl, aes(intercept = cc_sma_gspl[1,1], slope = cc_sma_gspl[2,1])) +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_gspl[2,1], digits = 2), sep = " "), size = 3, col = "black") +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_gspl[1,1], digits = 2), sep = " "), size = 3, col = "black")

#plot with outlier shown
plot_gspl_out <- plot_gspl +
  geom_point(aes(y=round(1.71,1), x=round(21.98,1)), colour="red", pch = 18) +
  ggtitle("Gardinerosergia splendens + outlier")
  
#plot
ggplotly(plot_gspl)
#plot with outlier
ggplotly(plot_gspl_out)

Neosergestes edwardsii

Note: this species has a big outlier that was removed for analyses and plots (eye 0.9, body 0.7). Also note that for this species, sample size was insufficient to estimate eye-body allometry (N = 6, p = 0.57 for correlation between eye size and body size).

# Subset data -----
nedw <- specimens %>% 
  filter(genus_species == "Neosergestes_edwardsii") %>%
  filter(!(Eye_Diameter == round(0.85,1) & Body_Length == round(0.72,1)))

# Fit SMA model ------
sma_nedw <- sma(formula = Eye_Diameter ~ Body_Length, 
            data = nedw, 
            log = "xy", #sets both x and y variables as logged
            method="SMA", #defines SMA as model
            alpha = 0.05)

#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_nedw, which = "default",type = "o") 
plot(sma_nedw, which = "residual",type = "o")
plot(sma_nedw, which = "qq", type = "o")
par(mfrow = c(1,1)) 

#view fits
summary(sma_nedw)
## Call: sma(formula = Eye_Diameter ~ Body_Length, data = nedw, log = "xy", 
##     method = "SMA", alpha = 0.05) 
## 
## Fit using Standardized Major Axis 
## 
## These variables were log-transformed before fitting: xy 
## 
## Confidence intervals (CI) are at 95%
## 
## ------------------------------------------------------------
## Coefficients:
##              elevation      slope
## estimate     1.4244363 -1.2712830
## lower limit -0.7533312 -3.8007834
## upper limit  3.6022037 -0.4252177
## 
## H0 : variables uncorrelated
## R-squared : 0.08540156 
## P-value : 0.57413
#save coefficients of fits as object
cc_sma_nedw <- data.frame(coef(sma_nedw))

# Make plot ------
plot_nedw <- ggplot(nedw, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) + 
  geom_point(size = 2, alpha = 0.8) + 
  scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
  scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
  scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
  ggtitle("Neosergestes edwardsii") + 
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) 

#plot with outlier shown
plot_nedw_out <- plot_nedw +
  geom_point(aes(x=round(0.85,1), y=round(0.72,1)), colour="red", pch = 18) +
  ggtitle("Neosergestes edwardsii + outlier")
  
#plot
ggplotly(plot_nedw)
#plot with outlier
ggplotly(plot_nedw_out)

Parasergestes armatus

# Subset data -----
parm <- specimens %>% filter(genus_species == "Parasergestes_armatus")

# Fit SMA model ------
sma_parm <- sma(formula = Eye_Diameter ~ Body_Length, 
            data = parm, 
            log = "xy", #sets both x and y variables as logged
            method="SMA", #defines SMA as model
            alpha = 0.05)

#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_parm, which = "default",type = "o") 
plot(sma_parm, which = "residual",type = "o")
plot(sma_parm, which = "qq", type = "o")
par(mfrow = c(1,1)) 

#view fits
summary(sma_parm)
## Call: sma(formula = Eye_Diameter ~ Body_Length, data = parm, log = "xy", 
##     method = "SMA", alpha = 0.05) 
## 
## Fit using Standardized Major Axis 
## 
## These variables were log-transformed before fitting: xy 
## 
## Confidence intervals (CI) are at 95%
## 
## ------------------------------------------------------------
## Coefficients:
##             elevation     slope
## estimate    -1.697577 1.0573117
## lower limit -2.095490 0.8217616
## upper limit -1.299663 1.3603801
## 
## H0 : variables uncorrelated
## R-squared : 0.5670657 
## P-value : 1.5744e-06
#save coefficients of fits as object
cc_sma_parm <- data.frame(coef(sma_parm))

# Make plot ------
plot_parm <- ggplot(parm, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) + 
  geom_point(size = 2, alpha = 0.8) + 
  scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
  scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
  scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
  ggtitle("Parasergestes armatus") + 
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
  geom_abline(data = cc_sma_parm, aes(intercept = cc_sma_parm[1,1], slope = cc_sma_parm[2,1])) +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_parm[2,1], digits = 2), sep = " "), size = 3, col = "black") +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_parm[1,1], digits = 2), sep = " "), size = 3, col = "black")

ggplotly(plot_parm)

Parasergestes vigilax

# Subset data -----
pvig <- specimens %>% filter(genus_species == "Parasergestes_vigilax")

# Fit SMA model ------
sma_pvig <- sma(formula = Eye_Diameter ~ Body_Length, 
            data = pvig, 
            log = "xy", #sets both x and y variables as logged
            method="SMA", #defines SMA as model
            alpha = 0.05)

#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_pvig, which = "default",type = "o") 
plot(sma_pvig, which = "residual",type = "o")
plot(sma_pvig, which = "qq", type = "o")
par(mfrow = c(1,1)) 

#view fits
summary(sma_pvig)
## Call: sma(formula = Eye_Diameter ~ Body_Length, data = pvig, log = "xy", 
##     method = "SMA", alpha = 0.05) 
## 
## Fit using Standardized Major Axis 
## 
## These variables were log-transformed before fitting: xy 
## 
## Confidence intervals (CI) are at 95%
## 
## ------------------------------------------------------------
## Coefficients:
##             elevation    slope
## estimate    -3.306489 2.375713
## lower limit -4.455793 1.641298
## upper limit -2.157186 3.438750
## 
## H0 : variables uncorrelated
## R-squared : 0.5274941 
## P-value : 0.00096145
#save coefficients of fits as object
cc_sma_pvig <- data.frame(coef(sma_pvig))

# Make plot ------
plot_pvig <- ggplot(pvig, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) + 
  geom_point(size = 2, alpha = 0.8) + 
  scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
  scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
  scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
  ggtitle("Parasergestes vigilax") + 
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
  geom_abline(data = cc_sma_pvig, aes(intercept = cc_sma_pvig[1,1], slope = cc_sma_pvig[2,1])) +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_pvig[2,1], digits = 2), sep = " "), size = 3, col = "black") +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_pvig[1,1], digits = 2), sep = " "), size = 3, col = "black")

ggplotly(plot_pvig)

Phorcosergia grandis

# Subset data -----
pgra <- specimens %>% filter(genus_species == "Phorcosergia_grandis")

# Fit SMA model ------
sma_pgra <- sma(formula = Eye_Diameter ~ Body_Length, 
            data = pgra, 
            log = "xy", #sets both x and y variables as logged
            method="SMA", #defines SMA as model
            alpha = 0.05)

#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_pgra, which = "default",type = "o") 
plot(sma_pgra, which = "residual",type = "o")
plot(sma_pgra, which = "qq", type = "o")
par(mfrow = c(1,1)) 

#view fits
summary(sma_pgra)
## Call: sma(formula = Eye_Diameter ~ Body_Length, data = pgra, log = "xy", 
##     method = "SMA", alpha = 0.05) 
## 
## Fit using Standardized Major Axis 
## 
## These variables were log-transformed before fitting: xy 
## 
## Confidence intervals (CI) are at 95%
## 
## ------------------------------------------------------------
## Coefficients:
##             elevation     slope
## estimate    -1.470665 0.9366565
## lower limit -1.632022 0.8480626
## upper limit -1.309308 1.0345055
## 
## H0 : variables uncorrelated
## R-squared : 0.8798127 
## P-value : < 2.22e-16
#save coefficients of fits as object
cc_sma_pgra <- data.frame(coef(sma_pgra))

# Make plot ------
plot_pgra <- ggplot(pgra, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) + 
  geom_point(size = 2, alpha = 0.8) + 
  scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
  scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
  scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
  ggtitle("Phorcosergia grandis") + 
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
  geom_abline(data = cc_sma_pgra, aes(intercept = cc_sma_pgra[1,1], slope = cc_sma_pgra[2,1])) +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_pgra[2,1], digits = 2), sep = " "), size = 3, col = "black") +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_pgra[1,1], digits = 2), sep = " "), size = 3, col = "black")

ggplotly(plot_pgra)

Robustasergia robusta

# Subset data -----
rrob <- specimens %>% filter(genus_species == "Robustasergia_robusta")

# Fit SMA model ------
sma_rrob <- sma(formula = Eye_Diameter ~ Body_Length, 
            data = rrob, 
            log = "xy", #sets both x and y variables as logged
            method="SMA", #defines SMA as model
            alpha = 0.05)

#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_rrob, which = "default",type = "o") 
plot(sma_rrob, which = "residual",type = "o")
plot(sma_rrob, which = "qq", type = "o")
par(mfrow = c(1,1)) 

#view fits
summary(sma_rrob)
## Call: sma(formula = Eye_Diameter ~ Body_Length, data = rrob, log = "xy", 
##     method = "SMA", alpha = 0.05) 
## 
## Fit using Standardized Major Axis 
## 
## These variables were log-transformed before fitting: xy 
## 
## Confidence intervals (CI) are at 95%
## 
## ------------------------------------------------------------
## Coefficients:
##              elevation     slope
## estimate    -1.2073093 0.8386026
## lower limit -1.5732188 0.6553974
## upper limit -0.8413998 1.0730198
## 
## H0 : variables uncorrelated
## R-squared : 0.6989428 
## P-value : 6.7833e-07
#save coefficients of fits as object
cc_sma_rrob <- data.frame(coef(sma_rrob))

# Make plot ------
plot_rrob <- ggplot(rrob, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) + 
  geom_point(size = 2, alpha = 0.8) + 
  scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
  scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
  scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
  ggtitle("Robustasergia robusta") + 
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
  geom_abline(data = cc_sma_rrob, aes(intercept = cc_sma_rrob[1,1], slope = cc_sma_rrob[2,1])) +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_rrob[2,1], digits = 2), sep = " "), size = 3, col = "black") +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_rrob[1,1], digits = 2), sep = " "), size = 3, col = "black")

ggplotly(plot_rrob)

Robustosergia regalis

# Subset data -----
rreg <- specimens %>% filter(genus_species == "Robustosergia_regalis")

# Fit SMA model ------
sma_rreg <- sma(formula = Eye_Diameter ~ Body_Length, 
            data = rreg, 
            log = "xy", #sets both x and y variables as logged
            method="SMA", #defines SMA as model
            alpha = 0.05)

#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_rreg, which = "default",type = "o") 
plot(sma_rreg, which = "residual",type = "o")
plot(sma_rreg, which = "qq", type = "o")
par(mfrow = c(1,1)) 

#view fits
summary(sma_rreg)
## Call: sma(formula = Eye_Diameter ~ Body_Length, data = rreg, log = "xy", 
##     method = "SMA", alpha = 0.05) 
## 
## Fit using Standardized Major Axis 
## 
## These variables were log-transformed before fitting: xy 
## 
## Confidence intervals (CI) are at 95%
## 
## ------------------------------------------------------------
## Coefficients:
##             elevation     slope
## estimate    -1.324440 0.8978007
## lower limit -1.531510 0.7807317
## upper limit -1.117369 1.0324239
## 
## H0 : variables uncorrelated
## R-squared : 0.8280327 
## P-value : 2.4971e-15
#save coefficients of fits as object
cc_sma_rreg <- data.frame(coef(sma_rreg))

# Make plot ------
plot_rreg <- ggplot(rreg, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) + 
  geom_point(size = 2, alpha = 0.8) + 
  scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
  scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
  scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
  ggtitle("Robustosergia regalis") + 
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
  geom_abline(data = cc_sma_rreg, aes(intercept = cc_sma_rreg[1,1], slope = cc_sma_rreg[2,1])) +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_rreg[2,1], digits = 2), sep = " "), size = 3, col = "black") +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_rreg[1,1], digits = 2), sep = " "), size = 3, col = "black")

ggplotly(plot_rreg)

Sergestes atlanticus

# Subset data -----
satl <- specimens %>% filter(genus_species == "Sergestes_atlanticus")

# Fit SMA model ------
sma_satl <- sma(formula = Eye_Diameter ~ Body_Length, 
            data = satl, 
            log = "xy", #sets both x and y variables as logged
            method="SMA", #defines SMA as model
            alpha = 0.05)

#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_satl, which = "default",type = "o") 
plot(sma_satl, which = "residual",type = "o")
plot(sma_satl, which = "qq", type = "o")
par(mfrow = c(1,1)) 

#view fits
summary(sma_satl)
## Call: sma(formula = Eye_Diameter ~ Body_Length, data = satl, log = "xy", 
##     method = "SMA", alpha = 0.05) 
## 
## Fit using Standardized Major Axis 
## 
## These variables were log-transformed before fitting: xy 
## 
## Confidence intervals (CI) are at 95%
## 
## ------------------------------------------------------------
## Coefficients:
##             elevation    slope
## estimate    -2.839127 1.951572
## lower limit -3.727822 1.418514
## upper limit -1.950431 2.684945
## 
## H0 : variables uncorrelated
## R-squared : 0.7067609 
## P-value : 8.6609e-05
#save coefficients of fits as object
cc_sma_satl <- data.frame(coef(sma_satl))

# Make plot ------
plot_satl <- ggplot(satl, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) + 
  geom_point(size = 2, alpha = 0.8) + 
  scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
  scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
  scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
  ggtitle("Sergestes atlanticus") + 
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
  geom_abline(data = cc_sma_satl, aes(intercept = cc_sma_satl[1,1], slope = cc_sma_satl[2,1])) +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_satl[2,1], digits = 2), sep = " "), size = 3, col = "black") +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_satl[1,1], digits = 2), sep = " "), size = 3, col = "black")

ggplotly(plot_satl)

Sergia tenuiremis

# Subset data -----
sten <- specimens %>% filter(genus_species == "Sergia_tenuiremis")

# Fit SMA model ------
sma_sten <- sma(formula = Eye_Diameter ~ Body_Length, 
            data = sten, 
            log = "xy", #sets both x and y variables as logged
            method="SMA", #defines SMA as model
            alpha = 0.05)

#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_sten, which = "default",type = "o") 
plot(sma_sten, which = "residual",type = "o")
plot(sma_sten, which = "qq", type = "o")
par(mfrow = c(1,1)) 

#view fits
summary(sma_sten)
## Call: sma(formula = Eye_Diameter ~ Body_Length, data = sten, log = "xy", 
##     method = "SMA", alpha = 0.05) 
## 
## Fit using Standardized Major Axis 
## 
## These variables were log-transformed before fitting: xy 
## 
## Confidence intervals (CI) are at 95%
## 
## ------------------------------------------------------------
## Coefficients:
##              elevation     slope
## estimate    -1.2313110 0.8025570
## lower limit -1.6822103 0.5747311
## upper limit -0.7804117 1.1206941
## 
## H0 : variables uncorrelated
## R-squared : 0.2880646 
## P-value : 0.0032339
#save coefficients of fits as object
cc_sma_sten <- data.frame(coef(sma_sten))

# Make plot ------
plot_sten <- ggplot(sten, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) + 
  geom_point(size = 2, alpha = 0.8) + 
  scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
  scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
  scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
  ggtitle("Sergia tenuiremis") + 
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
  geom_abline(data = cc_sma_sten, aes(intercept = cc_sma_sten[1,1], slope = cc_sma_sten[2,1])) +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_sten[2,1], digits = 2), sep = " "), size = 3, col = "black") +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_sten[1,1], digits = 2), sep = " "), size = 3, col = "black")

ggplotly(plot_sten)

Figure for static eye-body scaling across sergestid species

#make figure panels
fig.a <- plot_dcor + theme(legend.position = "none", axis.title.x = element_blank())
fig.b <- plot_dhen + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank())
fig.c <- plot_apec + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank())
fig.d <- plot_asar + theme(legend.position = "none", axis.title.x = element_blank())
fig.e <- plot_satl + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank())
fig.f <- plot_pvig + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank())
fig.g <- plot_parm + theme(legend.position = "none", axis.title.x = element_blank())
fig.h <- plot_euar + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank())
fig.i <- plot_gspl + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank())
fig.j <- plot_rreg + theme(legend.position = "none", axis.title.x = element_blank())
fig.k <- plot_rrob + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank())
fig.l <- plot_pgra + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank())
fig.m <- plot_sten + theme(legend.position = "none")
fig.n <- plot_ctal  + theme(legend.position = "none", axis.title.y = element_blank())
fig.o <- plot_chan + theme(legend.position = "none", axis.title.y = element_blank())

# arrange panels in figure
plots <- plot_grid(fig.a, fig.b, fig.c, fig.d, fig.e, fig.f, fig.g, fig.h, fig.i, fig.j, fig.k, fig.l, fig.m, fig.n, fig.o, #list of plots to arrange in grid
           align = 'vh', #align horizontally and vertically
           axis = 'lb',
           labels = c("A", "B", "C", "D", "E", "F", "G", "H", "I","J", "K", "L", "M", "N", "O"), #panel labels for figure
           hjust = -1, #adjustment for panel labels
           nrow = 5) #number of rows in grids

#extract legend
leg <- get_legend(fig.a + theme(legend.position="bottom"))


#view figure
plot_grid(plots, leg, nrow = 2, rel_heights = c(1, .01), align = 'vh')

#export figure
#png("../Figures/dev-allometry.png", width = 10, height = 12, units = "in", res = 500)
pdf("../Figures/static-allometry.pdf", width = 10, height = 12)
plot_grid(plots, leg, nrow = 2, rel_heights = c(1, .05), align = 'vh')
dev.off()

Test whether species differ in static allometry

We’ve looked at how eyes scale with body size within species, but are there significant differences in static allometry across species? Here we use two approaches: first, we compare linear mixed models with intercepts varying between species and slopes (+ intercepts) varying between species to see what best describes our data. Then, we look at pairwise comparisons among species to see which are driving these differences.

Linear mixed model to test whether slopes vary across species

Here, we use linear mixed models in the lme4 v.1.1.25 package to test species show differences in static eye-body allometric slopes. We fit 1) a constant-slope, variable intercept model to our eye and body size data with species as a random effect and 2) a variable slope model with the same effects, and then compared model fits by AIC scores. Both models are fit using reduced maximum likelihood (REML), as this approach is better than maximum likelihood (ML) for detecting differences in random effects among models, which is what we are testing here (ML is better for fixed effects).

First we removed the three outliers and the species with insufficient sampling (N. edwardsii) from our full specimen dataset, then ran each model.

#Remove outliers from specimen dataset
specimens_ed <- specimens %>%
  filter(!(genus_species == "Allosergestes_sargassi" & Eye_Diameter == 0.7 & Body_Length == 13.0)) %>% #(0.67, 12.97 raw)
  filter(!(genus_species == "Gardinerosergia_splendens" & Eye_Diameter == 1.7 & Body_Length == 22.0)) %>% #(1.71, 21.98 raw)
  filter(!(genus_species == "Neosergestes_edwardsii" & Eye_Diameter == 0.8 & Body_Length == 0.7)) #(0.85, 0.72 raw)

#Remove species with insufficient sampling
specimens_test <- specimens_ed %>%
  filter(genus_species != "Neosergestes_edwardsii") %>%
  mutate(species_text = as.factor(paste(Genus, Species, sep = " "))) # Add labels for species text

fixed slope model

# variable intercepts model (species can have different intercepts but share mean slope) -------
lme_fixed <- lmer(log10(Eye_Diameter) ~ log10(Body_Length) + (1|genus_species), 
                  data = specimens_test)

# summary 
summary(lme_fixed)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log10(Eye_Diameter) ~ log10(Body_Length) + (1 | genus_species)
##    Data: specimens_test
## 
## REML criterion at convergence: -1091.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.6547 -0.5660  0.0425  0.5699  3.7587 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  genus_species (Intercept) 0.004257 0.06525 
##  Residual                  0.004549 0.06745 
## Number of obs: 450, groups:  genus_species, 15
## 
## Fixed effects:
##                    Estimate Std. Error t value
## (Intercept)        -1.24683    0.05121  -24.35
## log10(Body_Length)  0.80642    0.03122   25.83
## 
## Correlation of Fixed Effects:
##             (Intr)
## lg10(Bdy_L) -0.941
# intercepts (variable) and slope (constant)
coef(lme_fixed)
## $genus_species
##                              (Intercept) log10(Body_Length)
## Allosergestes_pectinatus       -1.337988          0.8064169
## Allosergestes_sargassi         -1.295962          0.8064169
## Challengerosergia_hansjacobi   -1.164606          0.8064169
## Challengerosergia_talismani    -1.216875          0.8064169
## Deosergestes_corniculum        -1.328655          0.8064169
## Deosergestes_henseni           -1.278582          0.8064169
## Eusergestes_arcticus           -1.246866          0.8064169
## Gardinerosergia_splendens      -1.158266          0.8064169
## Parasergestes_armatus          -1.325352          0.8064169
## Parasergestes_vigilax          -1.299386          0.8064169
## Phorcosergia_grandis           -1.246064          0.8064169
## Robustasergia_robusta          -1.155267          0.8064169
## Robustosergia_regalis          -1.176361          0.8064169
## Sergestes_atlanticus           -1.234190          0.8064169
## Sergia_tenuiremis              -1.238014          0.8064169
## 
## attr(,"class")
## [1] "coef.mer"
# AIC
AIC(lme_fixed)
## [1] -1083.888
#fixed effects
fixef(lme_fixed)
##        (Intercept) log10(Body_Length) 
##         -1.2468288          0.8064169
#confidence intervals
## if these overlap zero, we do not have good support for significant effects
## if .sig01, which is our estimate of the variability in the random effect is very large and very widely defined, t indicates we may have a lack of precision between our groups - either because the group effect is small between groups, we have too few groups to get a more precise estimate, we have too few units within each group, or a combination of all of the above.
confint(lme_fixed, level = 0.95)
##                          2.5 %      97.5 %
## .sig01              0.04428477  0.09578174
## .sigma              0.06313279  0.07211335
## (Intercept)        -1.34880013 -1.14643360
## log10(Body_Length)  0.74536494  0.86929119

variable slope model

# variable slopes model (species can have different slopes/interaction effect between body size and species on eye size) -------
lme_variable <- lmer(log10(Eye_Diameter) ~ log10(Body_Length) + (1 + log10(Body_Length) | genus_species), 
                     data = specimens_test)

# summary
summary(lme_variable)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log10(Eye_Diameter) ~ log10(Body_Length) + (1 + log10(Body_Length) |  
##     genus_species)
##    Data: specimens_test
## 
## REML criterion at convergence: -1125.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.4279 -0.5359  0.0306  0.5910  4.2687 
## 
## Random effects:
##  Groups        Name               Variance Std.Dev. Corr 
##  genus_species (Intercept)        0.232807 0.48250       
##                log10(Body_Length) 0.104532 0.32331  -0.99
##  Residual                         0.003977 0.06306       
## Number of obs: 450, groups:  genus_species, 15
## 
## Fixed effects:
##                    Estimate Std. Error t value
## (Intercept)        -1.35693    0.13739  -9.876
## log10(Body_Length)  0.89778    0.09181   9.778
## 
## Correlation of Fixed Effects:
##             (Intr)
## lg10(Bdy_L) -0.993
# intercepts (variable) and slope (constant)
coef(lme_variable)
## $genus_species
##                              (Intercept) log10(Body_Length)
## Allosergestes_pectinatus      -2.2120964          1.4762026
## Allosergestes_sargassi        -1.7513868          1.1366814
## Challengerosergia_hansjacobi  -1.0760213          0.7484245
## Challengerosergia_talismani   -1.4111004          0.9318641
## Deosergestes_corniculum       -1.1953554          0.7273404
## Deosergestes_henseni          -0.7843247          0.4963146
## Eusergestes_arcticus          -1.2644263          0.8196651
## Gardinerosergia_splendens     -0.8716955          0.6169987
## Parasergestes_armatus         -1.3258605          0.8075505
## Parasergestes_vigilax         -2.0269469          1.3719045
## Phorcosergia_grandis          -1.3536721          0.8693323
## Robustasergia_robusta         -1.1498050          0.8037740
## Robustosergia_regalis         -1.2328226          0.8412590
## Sergestes_atlanticus          -1.8945693          1.2766884
## Sergia_tenuiremis             -0.8038697          0.5427274
## 
## attr(,"class")
## [1] "coef.mer"
# AIC
AIC(lme_variable)
## [1] -1113.879
#fixed effects
fixef(lme_variable)
##        (Intercept) log10(Body_Length) 
##         -1.3569302          0.8977818
#confidence intervals
## if these overlap zero, we do not have good support for significant effects
## if .sig01, which is our estimate of the variability in the random effect is very large and very widely defined, t indicates we may have a lack of precision between our groups - either because the group effect is small between groups, we have too few groups to get a more precise estimate, we have too few units within each group, or a combination of all of the above.
confint(lme_variable, level = 0.95)
##                          2.5 %      97.5 %
## .sig01              0.27582253  0.74974105
## .sig02             -0.99580190 -0.97348803
## .sig03              0.17701120  0.50987550
## .sigma              0.05903897  0.06767603
## (Intercept)        -1.64344107 -1.07983292
## log10(Body_Length)  0.71455335  1.09409513

Note that the parameter estimates for slopes and intercepts for species groups are a bit different than those we estimated via OLS. That’s because here, the model is looking across all the data for all the species, and finding a common mean slope/intercept as well as the group effects. Rather than assuming that the species we sampled are all possible species, fitting species as a random effect assumes that the species in our study are some random sampling of the possible species that exist, and that there is some true population mean among species. Here we do see differences across species, but the mixed model pulls the more extreme species toward the group mean (this is all based on variances, both within and between groups/species here). I think this is more appropriate for directly answering the question of whether species are following the same slope for static allometry or not, but it could be argued either way whether these slope and intercept estimates or the ones fitted individually to each species by SMA are more appropriate. If you are trying to predict eye size from body size, I would think that the individual species fit is more appropriate, beause that is a different question. So I think having both methods in the paper is reasonable (each addressing a specific, and different, question), but I am new to mixed models!

compare models

Then we can compare the two models via AIC and log likelihood.

#model AIC comparison
AIC(lme_fixed, lme_variable)
##              df       AIC
## lme_fixed     4 -1083.888
## lme_variable  6 -1113.879
#log likelihood comparison
anova(lme_fixed, lme_variable, refit=FALSE)
## Data: specimens_test
## Models:
## lme_fixed: log10(Eye_Diameter) ~ log10(Body_Length) + (1 | genus_species)
## lme_variable: log10(Eye_Diameter) ~ log10(Body_Length) + (1 + log10(Body_Length) | 
## lme_variable:     genus_species)
##              npar     AIC     BIC logLik deviance  Chisq Df Pr(>Chisq)    
## lme_fixed       4 -1083.9 -1067.5 545.94  -1091.9                         
## lme_variable    6 -1113.9 -1089.2 562.94  -1125.9 33.991  2  4.159e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Here we find by comparing AIC scores and log likelihoods that the model with variable slopes and intercepts is better supported than the one with fixed slopes and variable intercepts. We have evidence here that species differ significantly in the slopes of their static allometries rather than following a shared trajectory.

Test for significant effect of preservation as a covariate

Here, we add preservation as a fixed effect into our linear mixed models of eye size vs. body size with species as a random effect. We have already selected the variable slopes model, so here we use a maximum likelihood approach (ML) to fit models that 1) don’t include preservation as a covariate and 2) does include preservation as a covariate and then compare model fits. Note that when we chatted, I mentioned putting in preservation type as a random effect. However, I have read that random effects need to have more than 5 levels, and preservation type has only 3, so by default it needs to be fitted as a fixed effect here.

model without preservation type

# variable slopes model (species can have different slopes/interaction effect between body size and species on eye size) -------
lme_nopres <- lmer(log10(Eye_Diameter) ~ log10(Body_Length) + (1 + log10(Body_Length) | genus_species), 
                   data = specimens_test,
                   REML = FALSE) #uses ML instead of REML

# summary
summary(lme_nopres)
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: log10(Eye_Diameter) ~ log10(Body_Length) + (1 + log10(Body_Length) |  
##     genus_species)
##    Data: specimens_test
## 
##      AIC      BIC   logLik deviance df.resid 
##  -1123.3  -1098.6    567.6  -1135.3      444 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.4474 -0.5397  0.0290  0.5871  4.2627 
## 
## Random effects:
##  Groups        Name               Variance Std.Dev. Corr 
##  genus_species (Intercept)        0.205896 0.45376       
##                log10(Body_Length) 0.091875 0.30311  -0.99
##  Residual                         0.003983 0.06311       
## Number of obs: 450, groups:  genus_species, 15
## 
## Fixed effects:
##                    Estimate Std. Error t value
## (Intercept)        -1.35193    0.13036  -10.37
## log10(Body_Length)  0.89371    0.08686   10.29
## 
## Correlation of Fixed Effects:
##             (Intr)
## lg10(Bdy_L) -0.993
# intercepts (variable) and slope (constant)
coef(lme_nopres)
## $genus_species
##                              (Intercept) log10(Body_Length)
## Allosergestes_pectinatus      -2.1804369          1.4518087
## Allosergestes_sargassi        -1.7432805          1.1308540
## Challengerosergia_hansjacobi  -1.0817261          0.7520710
## Challengerosergia_talismani   -1.4101468          0.9312295
## Deosergestes_corniculum       -1.1915659          0.7253709
## Deosergestes_henseni          -0.7877675          0.4985130
## Eusergestes_arcticus          -1.2643293          0.8197731
## Gardinerosergia_splendens     -0.8765113          0.6200664
## Parasergestes_armatus         -1.3286959          0.8097119
## Parasergestes_vigilax         -1.9939820          1.3460557
## Phorcosergia_grandis          -1.3527444          0.8688290
## Robustasergia_robusta         -1.1602353          0.8095809
## Robustosergia_regalis         -1.2356690          0.8429100
## Sergestes_atlanticus          -1.8546819          1.2481093
## Sergia_tenuiremis             -0.8172020          0.5507936
## 
## attr(,"class")
## [1] "coef.mer"
# AIC
AIC(lme_nopres)
## [1] -1123.278
#fixed effects
fixef(lme_nopres)
##        (Intercept) log10(Body_Length) 
##         -1.3519317          0.8937118
#confidence intervals
## if these overlap zero, we do not have good support for significant effects
## if .sig01, which is our estimate of the variability in the random effect is very large and very widely defined, t indicates we may have a lack of precision between our groups - either because the group effect is small between groups, we have too few groups to get a more precise estimate, we have too few units within each group, or a combination of all of the above.
confint(lme_nopres, level = 0.95)
##                          2.5 %      97.5 %
## .sig01              0.27582245  0.74974090
## .sig02             -0.99580205 -0.97348803
## .sig03              0.17701116  0.50987543
## .sigma              0.05903897  0.06767603
## (Intercept)        -1.64344108 -1.07983293
## log10(Body_Length)  0.71455335  1.09409513

model including preservation type

#reorder levels for preservation so that contrasts will be compared to fresh samples
specimens_test$Preservation <- factor(specimens_test$Preservation, 
                                      levels = c("fresh", "ethanol", "paraformaldehyde"))

# variable slopes model (species can have different slopes/interaction effect between body size and species on eye size) -------
lme_withpres <- lmer(log10(Eye_Diameter) ~ log10(Body_Length) + Preservation + (1 + log10(Body_Length) | genus_species), 
                   data = specimens_test,
                   REML = FALSE) #uses ML instead of REML

# summary
summary(lme_withpres)
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: log10(Eye_Diameter) ~ log10(Body_Length) + Preservation + (1 +  
##     log10(Body_Length) | genus_species)
##    Data: specimens_test
## 
##      AIC      BIC   logLik deviance df.resid 
##  -1152.9  -1120.0    584.5  -1168.9      442 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.7037 -0.5122  0.0295  0.5644  4.3676 
## 
## Random effects:
##  Groups        Name               Variance Std.Dev. Corr 
##  genus_species (Intercept)        0.167876 0.40973       
##                log10(Body_Length) 0.071632 0.26764  -0.99
##  Residual                         0.003701 0.06083       
## Number of obs: 450, groups:  genus_species, 15
## 
## Fixed effects:
##                               Estimate Std. Error t value
## (Intercept)                  -1.353719   0.119794 -11.300
## log10(Body_Length)            0.894613   0.078169  11.445
## Preservationethanol          -0.019607   0.008151  -2.405
## Preservationparaformaldehyde  0.045520   0.012167   3.741
## 
## Correlation of Fixed Effects:
##             (Intr) l10(B_ Prsrvtnt
## lg10(Bdy_L) -0.990                
## Prsrvtnthnl -0.057  0.022         
## Prsrvtnprfr -0.068  0.041  0.455
# AIC
AIC(lme_withpres)
## [1] -1152.918
#fixed effects
fixef(lme_withpres)
##                  (Intercept)           log10(Body_Length) 
##                  -1.35371863                   0.89461257 
##          Preservationethanol Preservationparaformaldehyde 
##                  -0.01960702                   0.04551984
#confidence intervals
## if these overlap zero, we do not have good support for significant effects
## if .sig01, which is our estimate of the variability in the random effect is very large and very widely defined, t indicates we may have a lack of precision between our groups - either because the group effect is small between groups, we have too few groups to get a more precise estimate, we have too few units within each group, or a combination of all of the above.
confint(lme_withpres, level = 0.95)
##                                    2.5 %      97.5 %
## .sig01                        0.25281798  0.67323346
## .sig02                       -0.99341603 -0.96585906
## .sig03                        0.15905142  0.44860338
## .sigma                        0.05692288  0.06521125
## (Intercept)                  -1.61782966 -1.10480110
## log10(Body_Length)            0.73401476  1.07261072
## Preservationethanol          -0.03573800 -0.00350345
## Preservationparaformaldehyde  0.02141721  0.06962976

compare models

Then we can compare the two models via AIC and log likelihood.

#model AIC comparison
AIC(lme_nopres, lme_withpres)
##              df       AIC
## lme_nopres    6 -1123.278
## lme_withpres  8 -1152.918
#log likelihood comparison
anova(lme_nopres, lme_withpres)
## Data: specimens_test
## Models:
## lme_nopres: log10(Eye_Diameter) ~ log10(Body_Length) + (1 + log10(Body_Length) | 
## lme_nopres:     genus_species)
## lme_withpres: log10(Eye_Diameter) ~ log10(Body_Length) + Preservation + (1 + 
## lme_withpres:     log10(Body_Length) | genus_species)
##              npar     AIC     BIC logLik deviance Chisq Df Pr(>Chisq)    
## lme_nopres      6 -1123.3 -1098.6 567.64  -1135.3                        
## lme_withpres    8 -1152.9 -1120.0 584.46  -1168.9 33.64  2  4.957e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Here we see that the model that includes preservation type as a covariate is significantly better supported. If you take a look at the parameter estimates from the model, you can see that ethanol preservation has a negative effect on intercept (eye size compared to body size), but paraformaldehyde has a positive effect. Note also that these effects are very small, so while they are significantly affecting our models, they seem to have minor effects. I would use this to argue why we don’t include preservation type in all our ecology ancovas; our power is already low due to small sample sizes and including preservation as a covariate would probalby overparamaterize our models.

Linear SMA regression for pairwise comparisons among species slopes

The smatr package also allows you to test whether groups exhibit significant differences in allometric scaling (slopes or intercepts). Below, we again use standardized major axis regression to model the scaling of eye diameter with body length across all specimens measured and test whether we see significant differences in ontogenetic scaling across species.

*Note: As we did not have sufficient sampling to estimate eye-body allometry for Neosergestes edwardsii, I’ve removed it from the specimen dataset before fitting this model.

This test shows that among 15 sergestid species (450 specimens), species significantly differ in their allometric scaling relationships between eye diameter and body length (p < 0.0001). It also yields a table of pairwise comparisons so you can see exactly which species differ significantly in slope and which do not.

# Model for static allometry across specimens with species and species*body size interaction as covariates
sma_species <- sma(formula = Eye_Diameter ~ Body_Length * genus_species, 
            data = specimens_test, 
            log = "xy", #sets both x and y variables as logged
            method="SMA", #defines SMA as model
            multcomp = TRUE, # adds pairwise comparisons between groups
            multcompmethod = "adjusted", # adjusts p-value for multiple comparisons
            alpha = 0.05) 

#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_species, which = "default",type = "o") 
plot(sma_species,which = "residual",type = "o")
plot(sma_species, which = "qq", type = "o")
par(mfrow = c(1,1)) 

#view fits
summary(sma_species)
## Call: sma(formula = Eye_Diameter ~ Body_Length * genus_species, data = specimens_test, 
##     log = "xy", method = "SMA", alpha = 0.05, multcomp = TRUE, 
##     multcompmethod = "adjusted") 
## 
## Fit using Standardized Major Axis 
## 
## These variables were log-transformed before fitting: xy 
## 
## Confidence intervals (CI) are at 95%
## 
## ------------------------------------------------------------
## Results of comparing lines among groups.
## 
## H0 : slopes are equal.
## Likelihood ratio statistic : 78.37 with 14 degrees of freedom
## P-value : 5.6738e-11 
## ------------------------------------------------------------
## 
## Results of multiple comparisons among groups.
## 
## Test for pair-wise difference in slope :
##                  genus_species_1              genus_species_2       Pval
## 1       Allosergestes_pectinatus       Allosergestes_sargassi 1.00000000
## 2       Allosergestes_pectinatus Challengerosergia_hansjacobi 0.10631965
## 3       Allosergestes_pectinatus  Challengerosergia_talismani 0.03771725
## 4       Allosergestes_pectinatus      Deosergestes_corniculum 0.00370375
## 5       Allosergestes_pectinatus         Deosergestes_henseni 0.00191135
## 6       Allosergestes_pectinatus         Eusergestes_arcticus 1.00000000
## 7       Allosergestes_pectinatus    Gardinerosergia_splendens 0.00417512
## 8       Allosergestes_pectinatus        Parasergestes_armatus 0.14662731
## 9       Allosergestes_pectinatus        Parasergestes_vigilax 1.00000000
## 10      Allosergestes_pectinatus         Phorcosergia_grandis 0.00484593
## 11      Allosergestes_pectinatus        Robustasergia_robusta 0.00386781
## 12      Allosergestes_pectinatus        Robustosergia_regalis 0.00330606
## 13      Allosergestes_pectinatus         Sergestes_atlanticus 1.00000000
## 14      Allosergestes_pectinatus            Sergia_tenuiremis 0.00977338
## 15        Allosergestes_sargassi Challengerosergia_hansjacobi 0.17642048
## 16        Allosergestes_sargassi  Challengerosergia_talismani 0.03323729
## 17        Allosergestes_sargassi      Deosergestes_corniculum 0.00219121
## 18        Allosergestes_sargassi         Deosergestes_henseni 0.00062265
## 19        Allosergestes_sargassi         Eusergestes_arcticus 1.00000000
## 20        Allosergestes_sargassi    Gardinerosergia_splendens 0.00279183
## 21        Allosergestes_sargassi        Parasergestes_armatus 0.25381226
## 22        Allosergestes_sargassi        Parasergestes_vigilax 1.00000000
## 23        Allosergestes_sargassi         Phorcosergia_grandis 0.00071873
## 24        Allosergestes_sargassi        Robustasergia_robusta 0.00340588
## 25        Allosergestes_sargassi        Robustosergia_regalis 0.00055799
## 26        Allosergestes_sargassi         Sergestes_atlanticus 1.00000000
## 27        Allosergestes_sargassi            Sergia_tenuiremis 0.01682182
## 28  Challengerosergia_hansjacobi  Challengerosergia_talismani 1.00000000
## 29  Challengerosergia_hansjacobi      Deosergestes_corniculum 1.00000000
## 30  Challengerosergia_hansjacobi         Deosergestes_henseni 0.99999989
## 31  Challengerosergia_hansjacobi         Eusergestes_arcticus 1.00000000
## 32  Challengerosergia_hansjacobi    Gardinerosergia_splendens 1.00000000
## 33  Challengerosergia_hansjacobi        Parasergestes_armatus 1.00000000
## 34  Challengerosergia_hansjacobi        Parasergestes_vigilax 0.04921452
## 35  Challengerosergia_hansjacobi         Phorcosergia_grandis 1.00000000
## 36  Challengerosergia_hansjacobi        Robustasergia_robusta 1.00000000
## 37  Challengerosergia_hansjacobi        Robustosergia_regalis 1.00000000
## 38  Challengerosergia_hansjacobi         Sergestes_atlanticus 0.25203409
## 39  Challengerosergia_hansjacobi            Sergia_tenuiremis 1.00000000
## 40   Challengerosergia_talismani      Deosergestes_corniculum 0.99999999
## 41   Challengerosergia_talismani         Deosergestes_henseni 0.99993757
## 42   Challengerosergia_talismani         Eusergestes_arcticus 1.00000000
## 43   Challengerosergia_talismani    Gardinerosergia_splendens 1.00000000
## 44   Challengerosergia_talismani        Parasergestes_armatus 1.00000000
## 45   Challengerosergia_talismani        Parasergestes_vigilax 0.02124209
## 46   Challengerosergia_talismani         Phorcosergia_grandis 1.00000000
## 47   Challengerosergia_talismani        Robustasergia_robusta 0.99999992
## 48   Challengerosergia_talismani        Robustosergia_regalis 1.00000000
## 49   Challengerosergia_talismani         Sergestes_atlanticus 0.11697344
## 50   Challengerosergia_talismani            Sergia_tenuiremis 0.99999999
## 51       Deosergestes_corniculum         Deosergestes_henseni 1.00000000
## 52       Deosergestes_corniculum         Eusergestes_arcticus 0.99999999
## 53       Deosergestes_corniculum    Gardinerosergia_splendens 1.00000000
## 54       Deosergestes_corniculum        Parasergestes_armatus 1.00000000
## 55       Deosergestes_corniculum        Parasergestes_vigilax 0.00236873
## 56       Deosergestes_corniculum         Phorcosergia_grandis 1.00000000
## 57       Deosergestes_corniculum        Robustasergia_robusta 1.00000000
## 58       Deosergestes_corniculum        Robustosergia_regalis 1.00000000
## 59       Deosergestes_corniculum         Sergestes_atlanticus 0.01350778
## 60       Deosergestes_corniculum            Sergia_tenuiremis 1.00000000
## 61          Deosergestes_henseni         Eusergestes_arcticus 0.99999959
## 62          Deosergestes_henseni    Gardinerosergia_splendens 1.00000000
## 63          Deosergestes_henseni        Parasergestes_armatus 0.99999978
## 64          Deosergestes_henseni        Parasergestes_vigilax 0.00128436
## 65          Deosergestes_henseni         Phorcosergia_grandis 1.00000000
## 66          Deosergestes_henseni        Robustasergia_robusta 1.00000000
## 67          Deosergestes_henseni        Robustosergia_regalis 1.00000000
## 68          Deosergestes_henseni         Sergestes_atlanticus 0.00731827
## 69          Deosergestes_henseni            Sergia_tenuiremis 1.00000000
## 70          Eusergestes_arcticus    Gardinerosergia_splendens 0.99999997
## 71          Eusergestes_arcticus        Parasergestes_armatus 1.00000000
## 72          Eusergestes_arcticus        Parasergestes_vigilax 0.99996770
## 73          Eusergestes_arcticus         Phorcosergia_grandis 1.00000000
## 74          Eusergestes_arcticus        Robustasergia_robusta 0.99999993
## 75          Eusergestes_arcticus        Robustosergia_regalis 1.00000000
## 76          Eusergestes_arcticus         Sergestes_atlanticus 1.00000000
## 77          Eusergestes_arcticus            Sergia_tenuiremis 0.99999911
## 78     Gardinerosergia_splendens        Parasergestes_armatus 1.00000000
## 79     Gardinerosergia_splendens        Parasergestes_vigilax 0.00266011
## 80     Gardinerosergia_splendens         Phorcosergia_grandis 1.00000000
## 81     Gardinerosergia_splendens        Robustasergia_robusta 1.00000000
## 82     Gardinerosergia_splendens        Robustosergia_regalis 1.00000000
## 83     Gardinerosergia_splendens         Sergestes_atlanticus 0.01499340
## 84     Gardinerosergia_splendens            Sergia_tenuiremis 1.00000000
## 85         Parasergestes_armatus        Parasergestes_vigilax 0.06386567
## 86         Parasergestes_armatus         Phorcosergia_grandis 1.00000000
## 87         Parasergestes_armatus        Robustasergia_robusta 1.00000000
## 88         Parasergestes_armatus        Robustosergia_regalis 1.00000000
## 89         Parasergestes_armatus         Sergestes_atlanticus 0.31517941
## 90         Parasergestes_armatus            Sergia_tenuiremis 1.00000000
## 91         Parasergestes_vigilax         Phorcosergia_grandis 0.00320670
## 92         Parasergestes_vigilax        Robustasergia_robusta 0.00240837
## 93         Parasergestes_vigilax        Robustosergia_regalis 0.00220257
## 94         Parasergestes_vigilax         Sergestes_atlanticus 1.00000000
## 95         Parasergestes_vigilax            Sergia_tenuiremis 0.00481337
## 96          Phorcosergia_grandis        Robustasergia_robusta 1.00000000
## 97          Phorcosergia_grandis        Robustosergia_regalis 1.00000000
## 98          Phorcosergia_grandis         Sergestes_atlanticus 0.01851043
## 99          Phorcosergia_grandis            Sergia_tenuiremis 1.00000000
## 100        Robustasergia_robusta        Robustosergia_regalis 1.00000000
## 101        Robustasergia_robusta         Sergestes_atlanticus 0.01358845
## 102        Robustasergia_robusta            Sergia_tenuiremis 1.00000000
## 103        Robustosergia_regalis         Sergestes_atlanticus 0.01271145
## 104        Robustosergia_regalis            Sergia_tenuiremis 1.00000000
## 105         Sergestes_atlanticus            Sergia_tenuiremis 0.02662028
##         TestStat
## 1   7.666232e-02
## 2   1.070239e+01
## 3   1.269776e+01
## 4   1.710663e+01
## 5   1.836679e+01
## 6   1.514543e+00
## 7   1.687873e+01
## 8   1.006766e+01
## 9   8.341304e-01
## 10  1.659541e+01
## 11  1.702417e+01
## 12  1.732280e+01
## 13  8.875795e-04
## 14  1.526295e+01
## 15  9.695878e+00
## 16  1.293861e+01
## 17  1.810626e+01
## 18  2.051006e+01
## 19  1.218474e+00
## 20  1.764469e+01
## 21  8.943571e+00
## 22  1.268181e+00
## 23  2.023537e+01
## 24  1.726619e+01
## 25  2.072006e+01
## 26  7.817618e-02
## 27  1.423248e+01
## 28  2.850713e-03
## 29  1.459201e+00
## 30  2.159382e+00
## 31  7.165782e-01
## 32  1.474497e+00
## 33  3.154991e-03
## 34  1.218978e+01
## 35  7.281210e-01
## 36  1.670749e+00
## 37  1.215204e+00
## 38  8.958462e+00
## 39  1.652888e+00
## 40  1.971056e+00
## 41  2.909088e+00
## 42  7.908634e-01
## 43  1.901819e+00
## 44  1.381453e-02
## 45  1.378961e+01
## 46  1.203054e+00
## 47  2.130775e+00
## 48  1.877144e+00
## 49  1.051505e+01
## 50  1.917202e+00
## 51  1.091006e-01
## 52  1.922892e+00
## 53  9.383239e-03
## 54  1.528670e+00
## 55  1.795778e+01
## 56  4.919054e-01
## 57  3.949879e-02
## 58  9.043887e-02
## 59  1.464885e+01
## 60  1.486450e-01
## 61  2.285326e+00
## 62  4.454856e-02
## 63  2.224689e+00
## 64  1.912552e+01
## 65  1.138002e+00
## 66  1.181780e-02
## 67  4.397057e-01
## 68  1.581216e+01
## 69  1.724888e-02
## 70  2.036867e+00
## 71  6.530716e-01
## 72  2.808168e+00
## 73  1.445758e+00
## 74  2.123385e+00
## 75  1.716984e+00
## 76  1.501074e+00
## 77  2.362523e+00
## 78  1.546238e+00
## 79  1.773675e+01
## 80  5.367445e-01
## 81  9.239058e-03
## 82  1.427310e-01
## 83  1.445085e+01
## 84  8.641470e-02
## 85  1.169009e+01
## 86  7.995621e-01
## 87  1.740929e+00
## 88  1.285411e+00
## 89  8.475706e+00
## 90  1.723105e+00
## 91  1.738089e+01
## 92  1.792616e+01
## 93  1.809640e+01
## 94  6.837798e-01
## 95  1.660823e+01
## 96  7.055402e-01
## 97  2.447817e-01
## 98  1.405093e+01
## 99  7.801888e-01
## 100 2.376709e-01
## 101 1.463755e+01
## 102 4.529580e-02
## 103 1.476415e+01
## 104 3.827565e-01
## 105 1.336087e+01
## 
## ------------------------------------------------------------
## Coefficients by group in variable "genus_species"
## 
## Group: Allosergestes_pectinatus 
##             elevation    slope
## estimate    -2.812398 1.940053
## lower limit -3.491769 1.485663
## upper limit -2.133027 2.533418
## 
## H0 : variables uncorrelated.
## R-squared : 0.8156794 
## P-value : 9.643e-06 
## 
## Group: Allosergestes_sargassi 
##             elevation    slope
## estimate    -2.728530 1.844597
## lower limit -3.402995 1.419494
## upper limit -2.054065 2.397008
## 
## H0 : variables uncorrelated.
## R-squared : 0.4037707 
## P-value : 2.3852e-05 
## 
## Group: Challengerosergia_hansjacobi 
##             elevation    slope
## estimate    -1.528299 1.047036
## lower limit -1.913161 0.824145
## upper limit -1.143436 1.330207
## 
## H0 : variables uncorrelated.
## R-squared : 0.4129122 
## P-value : 3.3885e-06 
## 
## Group: Challengerosergia_talismani 
##             elevation     slope
## estimate    -1.578063 1.0389917
## lower limit -1.840235 0.8845412
## upper limit -1.315891 1.2204109
## 
## H0 : variables uncorrelated.
## R-squared : 0.8120894 
## P-value : 2.0471e-12 
## 
## Group: Deosergestes_corniculum 
##             elevation     slope
## estimate    -1.432046 0.8653057
## lower limit -1.737583 0.7020666
## upper limit -1.126509 1.0664999
## 
## H0 : variables uncorrelated.
## R-squared : 0.8421211 
## P-value : 8.1768e-08 
## 
## Group: Deosergestes_henseni 
##             elevation     slope
## estimate    -1.306994 0.8238705
## lower limit -1.593941 0.6637531
## upper limit -1.020047 1.0226131
## 
## H0 : variables uncorrelated.
## R-squared : 0.3380975 
## P-value : 1.6982e-06 
## 
## Group: Eusergestes_arcticus 
##              elevation     slope
## estimate    -2.1539322 1.3306114
## lower limit -3.8635100 0.6694485
## upper limit -0.4443543 2.6447543
## 
## H0 : variables uncorrelated.
## R-squared : 0.8368056 
## P-value : 0.029484 
## 
## Group: Gardinerosergia_splendens 
##              elevation     slope
## estimate    -1.2253950 0.8523703
## lower limit -1.5326713 0.6732695
## upper limit -0.9181186 1.0791149
## 
## H0 : variables uncorrelated.
## R-squared : 0.4597308 
## P-value : 1.1132e-06 
## 
## Group: Parasergestes_armatus 
##             elevation     slope
## estimate    -1.697577 1.0573117
## lower limit -2.095490 0.8217616
## upper limit -1.299663 1.3603801
## 
## H0 : variables uncorrelated.
## R-squared : 0.5670657 
## P-value : 1.5744e-06 
## 
## Group: Parasergestes_vigilax 
##             elevation    slope
## estimate    -3.306489 2.375713
## lower limit -4.455793 1.641298
## upper limit -2.157186 3.438750
## 
## H0 : variables uncorrelated.
## R-squared : 0.5274941 
## P-value : 0.00096145 
## 
## Group: Phorcosergia_grandis 
##             elevation     slope
## estimate    -1.470665 0.9366565
## lower limit -1.632022 0.8480626
## upper limit -1.309308 1.0345055
## 
## H0 : variables uncorrelated.
## R-squared : 0.8798127 
## P-value : < 2.22e-16 
## 
## Group: Robustasergia_robusta 
##              elevation     slope
## estimate    -1.2073093 0.8386026
## lower limit -1.5732188 0.6553974
## upper limit -0.8413998 1.0730198
## 
## H0 : variables uncorrelated.
## R-squared : 0.6989428 
## P-value : 6.7833e-07 
## 
## Group: Robustosergia_regalis 
##             elevation     slope
## estimate    -1.324440 0.8978007
## lower limit -1.531510 0.7807317
## upper limit -1.117369 1.0324239
## 
## H0 : variables uncorrelated.
## R-squared : 0.8280327 
## P-value : 2.4971e-15 
## 
## Group: Sergestes_atlanticus 
##             elevation    slope
## estimate    -2.839127 1.951572
## lower limit -3.727822 1.418514
## upper limit -1.950431 2.684945
## 
## H0 : variables uncorrelated.
## R-squared : 0.7067609 
## P-value : 8.6609e-05 
## 
## Group: Sergia_tenuiremis 
##              elevation     slope
## estimate    -1.2313110 0.8025570
## lower limit -1.6822103 0.5747311
## upper limit -0.7804117 1.1206941
## 
## H0 : variables uncorrelated.
## R-squared : 0.2880646 
## P-value : 0.0032339
#save coefficients of fits as object
cc_species <- data.frame(coef(sma_species))

This figure shows the static allometry of all the species with sufficient sampling (n = 15). It’s interactive, so you can click on different species in the legend to hide or show those points. That can be helpful for exploring patterns in a big mess of data. To me, it looks a bit like the smaller-bodied species have higher slopes? It looks like Sergestes atlanticus, Allosergestes_pectinatus, and Parasergestes vigilax have different (higher) slopes compared to the other sergestids, and the other species differ a bit in slopes but more so in intercepts. You can probably pull out which species are driving the difference by looking across the pairwise comparisons.

#reorder factor levels for figure legend (phylo order)
specimens_test$species_text <- factor(specimens_test$species_text, 
                                      levels = c("Deosergestes corniculum",
                                                 "Deosergestes henseni",
                                                 "Allosergestes pectinatus",
                                                 "Allosergestes sargassi",
                                                 "Sergestes atlanticus",
                                                 "Parasergestes vigilax",
                                                 "Parasergestes armatus",
                                                 "Eusergestes arcticus",
                                                 "Gardinerosergia splendens",
                                                 "Robustosergia regalis",
                                                 "Robustasergia robusta",
                                                 "Phorcosergia grandis",
                                                 "Sergia tenuiremis",
                                                 "Challengerosergia talismani",
                                                 "Challengerosergia hansjacobi"))

#make shape pallette
shapes.sp <- c("Deosergestes corniculum" = 21,
              "Deosergestes henseni" = 22, 
              "Allosergestes pectinatus" = 23, 
              "Allosergestes sargassi" = 24,
              "Sergestes atlanticus" = 25, 
              "Parasergestes vigilax" = 21, 
              "Parasergestes armatus" = 22,
              "Eusergestes arcticus" = 23, 
              "Gardinerosergia splendens" = 24, 
              "Robustosergia regalis" = 25, 
              "Robustasergia robusta" = 21, 
              "Phorcosergia grandis" = 22,
              "Sergia tenuiremis" = 23,
              "Challengerosergia talismani" = 24,
              "Challengerosergia hansjacobi" = 25)

cols.sp <- c("Deosergestes corniculum" = "#16ABA8",
              "Deosergestes henseni" = "#1f78b4", 
              "Allosergestes pectinatus" = "#3acf75", 
              "Allosergestes sargassi" = "#33a02c",
              "Sergestes atlanticus" = "#fb9a99", 
              "Parasergestes vigilax" = "#e31a1c", 
              "Parasergestes armatus" = "#fdbf6f",
              "Eusergestes arcticus" = "#ff7f00", 
              "Gardinerosergia splendens" = "#cab2d6",
              "Robustosergia regalis" = "#6a3d9a", 
              "Robustasergia robusta" = "gray31",
              "Phorcosergia grandis" = "burlywood4",
              "Sergia tenuiremis" = "black",
              "Challengerosergia talismani" = "maroon1",
              "Challengerosergia hansjacobi" = "orchid1")


#set s limits by species for plots
limits <- specimens_test %>%
  group_by(genus_species) %>%
  summarise(max = max(Body_Length, na.rm = TRUE),
            min = min(Body_Length, na.rm = TRUE))

#make plot
plot_species <- ggplot(specimens_test, aes(x = Body_Length, y = Eye_Diameter, color =  species_text, shape = species_text, fill = species_text)) + 
  geom_point(size = 2, alpha = 0.35) + 
  scale_shape_manual(values = shapes.sp, name = "Species") + 
  scale_color_manual(values = cols.sp, name = "Species") +
  scale_fill_manual(values = cols.sp, name = "Species") +
  scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
  scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.text = element_text(face = "italic")) +
    geom_segment(aes(x = limits$min[which(limits$genus_species == "Allosergestes_pectinatus")], xend = limits$max[which(limits$genus_species == "Allosergestes_pectinatus")], y = 10^(log10(limits$min[which(limits$genus_species == "Allosergestes_pectinatus")])*cc_species["Allosergestes_pectinatus", "slope"] + cc_species["Allosergestes_pectinatus", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Allosergestes_pectinatus")])*cc_species["Allosergestes_pectinatus", "slope"] + cc_species["Allosergestes_pectinatus", "elevation"])), color = cols.sp["Allosergestes pectinatus"]) + #Rana temporaria adult scaling
  geom_segment(aes(x = limits$min[which(limits$genus_species == "Allosergestes_sargassi")], xend = limits$max[which(limits$genus_species == "Allosergestes_sargassi")], y = 10^(log10(limits$min[which(limits$genus_species == "Allosergestes_sargassi")])*cc_species["Allosergestes_sargassi", "slope"] + cc_species["Allosergestes_sargassi", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Allosergestes_sargassi")])*cc_species["Allosergestes_sargassi", "slope"] + cc_species["Allosergestes_sargassi", "elevation"])), color = cols.sp["Allosergestes sargassi"]) + #Allosergestes_sargassi adult scaling
  geom_segment(aes(x = limits$min[which(limits$genus_species == "Challengerosergia_hansjacobi")], xend = limits$max[which(limits$genus_species == "Challengerosergia_hansjacobi")], y = 10^(log10(limits$min[which(limits$genus_species == "Challengerosergia_hansjacobi")])*cc_species["Challengerosergia_hansjacobi", "slope"] + cc_species["Challengerosergia_hansjacobi", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Challengerosergia_hansjacobi")])*cc_species["Challengerosergia_hansjacobi", "slope"] + cc_species["Challengerosergia_hansjacobi", "elevation"])), color = cols.sp["Challengerosergia hansjacobi"]) + #Challengerosergia_hansjacobi adult scaling
  geom_segment(aes(x = limits$min[which(limits$genus_species == "Challengerosergia_talismani")], xend = limits$max[which(limits$genus_species == "Challengerosergia_talismani")], y = 10^(log10(limits$min[which(limits$genus_species == "Challengerosergia_talismani")])*cc_species["Challengerosergia_talismani", "slope"] + cc_species["Challengerosergia_talismani", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Challengerosergia_talismani")])*cc_species["Challengerosergia_talismani", "slope"] + cc_species["Challengerosergia_talismani", "elevation"])), color = cols.sp["Challengerosergia talismani"]) + #Challengerosergia_talismani adult scaling
  geom_segment(aes(x = limits$min[which(limits$genus_species == "Deosergestes_corniculum")], xend = limits$max[which(limits$genus_species == "Deosergestes_corniculum")], y = 10^(log10(limits$min[which(limits$genus_species == "Deosergestes_corniculum")])*cc_species["Deosergestes_corniculum", "slope"] + cc_species["Deosergestes_corniculum", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Deosergestes_corniculum")])*cc_species["Deosergestes_corniculum", "slope"] + cc_species["Deosergestes_corniculum", "elevation"])), color = cols.sp["Deosergestes corniculum"]) + #Deosergestes_corniculum adult scaling
  geom_segment(aes(x = limits$min[which(limits$genus_species == "Deosergestes_henseni")], xend = limits$max[which(limits$genus_species == "Deosergestes_henseni")], y = 10^(log10(limits$min[which(limits$genus_species == "Deosergestes_henseni")])*cc_species["Deosergestes_henseni", "slope"] + cc_species["Deosergestes_henseni", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Deosergestes_henseni")])*cc_species["Deosergestes_henseni", "slope"] + cc_species["Deosergestes_henseni", "elevation"])), color = cols.sp["Deosergestes henseni"]) + #Deosergestes henseni adult scaling
  geom_segment(aes(x = limits$min[which(limits$genus_species == "Gardinerosergia_splendens")], xend = limits$max[which(limits$genus_species == "Gardinerosergia_splendens")], y = 10^(log10(limits$min[which(limits$genus_species == "Gardinerosergia_splendens")])*cc_species["Gardinerosergia_splendens", "slope"] + cc_species["Gardinerosergia_splendens", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Gardinerosergia_splendens")])*cc_species["Gardinerosergia_splendens", "slope"] + cc_species["Gardinerosergia_splendens", "elevation"])), color = cols.sp["Gardinerosergia splendens"]) + #Gardinerosergia splendens adult scaling
  geom_segment(aes(x = limits$min[which(limits$genus_species == "Parasergestes_armatus")], xend = limits$max[which(limits$genus_species == "Parasergestes_armatus")], y = 10^(log10(limits$min[which(limits$genus_species == "Parasergestes_armatus")])*cc_species["Parasergestes_armatus", "slope"] + cc_species["Parasergestes_armatus", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Parasergestes_armatus")])*cc_species["Parasergestes_armatus", "slope"] + cc_species["Parasergestes_armatus", "elevation"])), color = cols.sp["Parasergestes armatus"]) + #Parasergestes armatus adult scaling
  geom_segment(aes(x = limits$min[which(limits$genus_species == "Parasergestes_vigilax")], xend = limits$max[which(limits$genus_species == "Parasergestes_vigilax")], y = 10^(log10(limits$min[which(limits$genus_species == "Parasergestes_vigilax")])*cc_species["Parasergestes_vigilax", "slope"] + cc_species["Parasergestes_vigilax", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Parasergestes_vigilax")])*cc_species["Parasergestes_vigilax", "slope"] + cc_species["Parasergestes_vigilax", "elevation"])), color = cols.sp["Parasergestes vigilax"]) + #Parasergestes vigilax adult scaling
  geom_segment(aes(x = limits$min[which(limits$genus_species == "Phorcosergia_grandis")], xend = limits$max[which(limits$genus_species == "Phorcosergia_grandis")], y = 10^(log10(limits$min[which(limits$genus_species == "Phorcosergia_grandis")])*cc_species["Phorcosergia_grandis", "slope"] + cc_species["Phorcosergia_grandis", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Phorcosergia_grandis")])*cc_species["Phorcosergia_grandis", "slope"] + cc_species["Phorcosergia_grandis", "elevation"])), color = cols.sp["Phorcosergia grandis"]) + #Phorcosergia grandis adult scaling
  geom_segment(aes(x = limits$min[which(limits$genus_species == "Robustasergia_robusta")], xend = limits$max[which(limits$genus_species == "Robustasergia_robusta")], y = 10^(log10(limits$min[which(limits$genus_species == "Robustasergia_robusta")])*cc_species["Robustasergia_robusta", "slope"] + cc_species["Robustasergia_robusta", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Robustasergia_robusta")])*cc_species["Robustasergia_robusta", "slope"] + cc_species["Robustasergia_robusta", "elevation"])), color = cols.sp["Robustasergia robusta"]) + #Robustasergia robusta adult scaling
  geom_segment(aes(x = limits$min[which(limits$genus_species == "Robustosergia_regalis")], xend = limits$max[which(limits$genus_species == "Robustosergia_regalis")], y = 10^(log10(limits$min[which(limits$genus_species == "Robustosergia_regalis")])*cc_species["Robustosergia_regalis", "slope"] + cc_species["Robustosergia_regalis", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Robustosergia_regalis")])*cc_species["Robustosergia_regalis", "slope"] + cc_species["Robustosergia_regalis", "elevation"])), color = cols.sp["Robustosergia regalis"]) +
  geom_segment(aes(x = limits$min[which(limits$genus_species == "Sergestes_atlanticus")], xend = limits$max[which(limits$genus_species == "Sergestes_atlanticus")], y = 10^(log10(limits$min[which(limits$genus_species == "Sergestes_atlanticus")])*cc_species["Sergestes_atlanticus", "slope"] + cc_species["Sergestes_atlanticus", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Sergestes_atlanticus")])*cc_species["Sergestes_atlanticus", "slope"] + cc_species["Sergestes_atlanticus", "elevation"])), color = cols.sp["Sergestes atlanticus"]) +
  geom_segment(aes(x = limits$min[which(limits$genus_species == "Sergia_tenuiremis")], xend = limits$max[which(limits$genus_species == "Sergia_tenuiremis")], y = 10^(log10(limits$min[which(limits$genus_species == "Sergia_tenuiremis")])*cc_species["Sergia_tenuiremis", "slope"] + cc_species["Sergia_tenuiremis", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Sergia_tenuiremis")])*cc_species["Sergia_tenuiremis", "slope"] + cc_species["Sergia_tenuiremis", "elevation"])), color = cols.sp["Sergia tenuiremis"])+
  geom_segment(aes(x = limits$min[which(limits$genus_species == "Eusergestes_arcticus")], xend = limits$max[which(limits$genus_species == "Eusergestes_arcticus")], y = 10^(log10(limits$min[which(limits$genus_species == "Eusergestes_arcticus")])*cc_species["Eusergestes_arcticus", "slope"] + cc_species["Eusergestes_arcticus", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Eusergestes_arcticus")])*cc_species["Eusergestes_arcticus", "slope"] + cc_species["Eusergestes_arcticus", "elevation"])), color = cols.sp["Eusergestes arcticus"])


#interactive plot
ggplotly(plot_species)
#make plot with alpha = 1 for legend
plot_species_leg <- ggplot(specimens_test, aes(x = Body_Length, y = Eye_Diameter, color =  species_text, shape = species_text, fill = species_text)) + 
  geom_point(size = 2, alpha = 1) + 
  scale_shape_manual(values = shapes.sp, name = "Species") + 
  scale_color_manual(values = cols.sp, name = "Species") +
  scale_fill_manual(values = cols.sp, name = "Species") +
  theme_bw() 

# remove legend from previous plot
plot <- plot_species + theme(legend.position = "none")

# extract legend from this plot with alpha = 1
leg <- get_legend(plot_species_leg)

# construct full figure with transparent fig and non-transparent legend
plot_grid(plot, leg, ncol = 2, rel_widths = c(1, .4))

#export figure
#png("../Figures/dev-species.png", width = 12, height = 12, units = "in", res = 500)
pdf("../Figures/static-species.pdf", width = 9, height = 5)
plot_grid(plot, leg, ncol = 2, rel_widths = c(1, .4))
dev.off()

Phylogenetic distribution of allometric slopes

Below, we plot species-specific slopes for developmental eye-body allometry onto the sergestid tree. There is potential to do more with these data; do we have predictions for what might drive species differences in developmental allometry? It looks like body size might be one factor, as the smaller-bodied species seem to have higher slopes. Whether species exhibit ontogenetic descent may be another factor to explore? We might predict that species that have ontogenetic descent and are going into dimmer habitats as they grow may invest in higher eye growth relative to body growth. Do all sergestids exhibit ontogenetic descent? So far, all I’ve done with these data in an evolutionary framework is test for phylogenetic signal. Let me know if you think of anything else cool to look at! Maybe we could test whether larger body sizes are correlated with lower slopes? Seemed from the plot of all the species that maybe the smaller bodied species had the higher slopes..

#Add slopes for developmental allometry into species dataframe ----

#add coefficient rownames to genus_species column
cc_species$genus_species <- rownames(cc_species)

#merge coefficients with species data by rowname
species_sma <- left_join(traits, cc_species, by = "genus_species") %>%
  filter(genus_species != "Neosergestes_edwardsii")

# Prep data and phylogeny -----

#Make list of taxa to drop (in tree but not in dataset)
drops <- setdiff(tree$tip.label, as.character(species_sma$Binomial))

#Drop unwanted tips from phylogeny
tree.dev <- drop.tip(phy = tree, tip = drops) 

# Make Genus_species tip labels to replace Binomial codes
sp.tips <- data.frame(old = species_sma$Binomial, new = species_sma$genus_species)

# replace tip labels in phylogeny
tree.dev.plot <- sub.taxa.label(tree.dev, sp.tips)

#change row names of species data
rownames(species_sma) <- species_sma$genus_species

#check that phylogeny and data match exactly
name.check(tree.dev.plot, species_sma)

#resort trait dataset to the order of tree tip labels
species_sma <- species_sma[tree.dev.plot$tip.label, ] 

#make trait vector for slope
slope <- as.vector(species_sma$slope) 

#add tip label names to vector
names(slope) <- species_sma$genus_species

# Phylogeny with slopes for developmental eye-body allometry----
plotTree.wBars(tree.dev.plot, slope, 
               scale = 0.5, 
               tip.labels = TRUE,
               offset = 0.1,
               ftype = "bi",
               fsize = 1)

#export figure
#png("../Figures/phylo-slopes.png", width = 10, height = 7, units = "in", res = 500)
pdf("../Figures/phylo-slopes.pdf", width = 10, height = 7)

plotTree.wBars(tree.dev.plot, slope, 
               scale = 0.5, 
               tip.labels = TRUE,
               offset = 0.1,
               ftype = "bi",
               fsize = 1)
dev.off()

To test for phylogenetic signal, we fit Pagel’s lambda using a maximum likelihood approach in caper v.1.0.1.

#fit Pagel's lambda in caper
slopes <- data.frame(slope, names(slope))
mod <- pgls(slope ~ 1, comparative.data(tree.dev.plot, slopes,"names.slope."), lambda="ML")

#summary of model
summary(mod)
## 
## Call:
## pgls(formula = slope ~ 1, data = comparative.data(tree.dev.plot, 
##     slopes, "names.slope."), lambda = "ML")
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.97761 -0.02774  0.04743  0.33125  0.87950 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 0.770
##    lower bound : 0.000, p = 0.20171
##    upper bound : 1.000, p = 0.63856
##    95.0% CI   : (NA, NA)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##             Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)  1.30332    0.19267  6.7644 9.105e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5198 on 14 degrees of freedom
## Multiple R-squared:     0,   Adjusted R-squared:     0 
## F-statistic:   NaN on 0 and 14 DF,  p-value: NA
#plot likelihood of lambda
plot(pgls.profile(mod, "lambda"))

The estimate for Pagel’s lambda is 0.77. However, this estimate was not significant, and neither the upper or lower bound of the confidence interval could be estimated. This could indicate that there is no phylogenetic signal in these data, but could also be the result of our sample size of species for comparative evolutionary models (n = 15). Simulations have previously found that estimating Pagel’s lambda for datasets of >20 species is problematic and often unreliable (see Freckleton et al., 2002). The likelihood profile of lambda shows that there is increased likelihood that it’s somewhere between 0.5 and 1.

Evolutionary eye-body allometry in Sergestidae

We next examine evolutionary (interspecific) eye-body allometry across sergestids using species means for eye diameter and body length. In calculating species means, we want to use adults only and also make sure that the data we use downstream for traits is relevant to the size class of individuals that we are taking species means from. The biggest concern for the Sergestidae is ontogenetic descent, as this is common across midwater animals and could confound relationships between eye size and depth. However, body size for ‘adult’ depth distribution data we used in this paper was cutoff between 4-11mm depending on species, and our smallest specimen here is 13 mm, so we are confident that species averages based on full species samples in our dataset are appropriate here.

Below, we use phylogenetic least-squares regression in caper to fit the allometric relationship between species means for eye size and body length across 16 species of sergestids.

# Prep data ------

#Calculate species means (3 outliers removed)
sp_means <- specimens_ed %>% 
  mutate_if(is.character, as.factor) %>% 
  group_by(genus_species) %>%
  summarise(eye_av = mean(Eye_Diameter), 
            length_av = mean(Body_Length), 
            n = n()) %>%
  ungroup()

#Merge with species trait data stored in dataframe 'species'
species <- traits %>%
  left_join(sp_means, by = "genus_species") %>%
  mutate(Organ = factor(Organ, levels = c("pesta","lensed", "unlensed", "none")))

#check that names match in dataframe and tree
name.check(phy = tree, data = species, data.names = species$Binomial)
  
#use caper function to combine phylogeny and data into one object (this function also matches species names in tree and dataset)
species.comp <- comparative.data(data = species, phy = tree, names.col = "Binomial", vcv = TRUE, na.omit = FALSE, warn.dropped = TRUE)

#check for dropped tips or dropped species
species.comp$dropped$tips #phylogeny
species.comp$dropped$unmatched.rows #dataset

# PGLS model -------

#fit model for eye diameter ~ body length
pgls_evo <- pgls(log10(eye_av) ~ log10(length_av), 
               data = species.comp, 
               lambda = "ML", #uses Maximum Liklihood estimate of lambda
               param.CI = 0.95)

#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_evo)

par(mfrow = c(1, 1))

The standard model checks for the PGLS of interspecific eye diameter vs. body length are a bit wonky. The residuals look a little bimodal when they should be normal (though don’t have huge weird outliers). We’ll have to decide if we’re ok with this. Assuming we accept this, model results are below.

#print model output 
summary(pgls_evo)
## 
## Call:
## pgls(formula = log10(eye_av) ~ log10(length_av), data = species.comp, 
##     lambda = "ML", param.CI = 0.95)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.100464 -0.041912  0.009141  0.044272  0.080878 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 1.000
##    lower bound : 0.000, p = 0.077147
##    upper bound : 1.000, p = 1    
##    95.0% CI   : (NA, NA)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##                   Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)      -1.288069   0.147956 -8.7057 5.051e-07 ***
## log10(length_av)  0.822852   0.095634  8.6042 5.810e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.05855 on 14 degrees of freedom
## Multiple R-squared: 0.841,   Adjusted R-squared: 0.8296 
## F-statistic: 74.03 on 1 and 14 DF,  p-value: 5.81e-07
#save coefficients of fits as object
cc_pgls_evo <- data.frame(coef(pgls_evo))

The PGLS regression shows that eyes scale with body length with a slope of 0.82, indicating that eyes have a negative evolutionary allometry (hypoallometric, slope < 1) with body length. This is what is most common among all vertebrates and some invertebrates that have been examined so far (inverts tend to vary a lot more than verts).

#Likelihood plot for Pagel's lambda 
lambda.evo <- pgls.profile(pgls_evo, "lambda")
plot(lambda.evo, main = "Pagel's lambda", 
     sub = "Solid red line indicates estimate for lambda and broken red lines indcaite the 95% confidence interval")

The maximum likelihood estimate of lambda is 1.00 (high phylogenetic signal/Brownian Motion), so that’s what caper has used to create the variance-covariance matrix for the regression. However, small sample sizes (< 20-30 data points) have no/little power to estimate lambda (see Freckleton et al. 2002), so tests for phylogenetic signal in this dataset are essentially meaningless (note that the 95% confidence interval for lambda here is 0 to 1, all possible values, and the p-values are non-signficant). For this reason, if any of our other comparative analyses estimate a lambda of 0, we should repeat them again forcing lambda to 1, just so we can see what difference that makes in the model fits (as we can’t know what lambda actually is in this small dataset).

However, if we choose to keep this model, we can plot the model fit onto our data.

# Make plot 
plot_evo <- ggplot(species, aes(x = Body_length, y = Eye_diameter, text = genus_species)) + 
  geom_point(size = 2, alpha = 0.6) + 
  scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
  scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  geom_abline(data = cc_pgls_evo, aes(intercept = cc_pgls_evo[1,1], slope = cc_pgls_evo[2,1])) +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_pgls_evo[2,1], digits = 2), sep = " "), size = 3) +
  geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_pgls_evo[1,1], digits = 2), sep = " "), size = 3)

# Interactive plot
ggplotly(plot_evo)

Bounded extremes of models

With sample sizes too small to generate a robust estimate of lambda, one option is to run models with the lowest and possible bounds of lambda (0 and 1), so that we can see what the range of possible parameter estimates are depending on the phylogenetic signal of the model residuals. Below, we fit these two models so that we can compare them.

Here is the model fit with lambda fixed at 0 (no phylogenetic signal in model residuals).

#fit model for eye diameter ~ body length (lambda = 0)
pgls_evo0 <- pgls(log10(eye_av) ~ log10(length_av), 
               data = species.comp, 
               lambda = 0, #set lambda to 0
               bounds=list(lambda=c(0,0)),
               param.CI = 0.95)

#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_evo0)

par(mfrow = c(1, 1))

#parameter estimates
summary(pgls_evo0)
## 
## Call:
## pgls(formula = log10(eye_av) ~ log10(length_av), data = species.comp, 
##     lambda = 0, param.CI = 0.95, bounds = list(lambda = c(0, 
##         0)))
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.108877 -0.031123 -0.008501  0.040904  0.086975 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [Fix]  : 0.000
## delta  [Fix]  : 1.000
## 
## Coefficients:
##                   Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)      -1.431224   0.141272 -10.131 7.949e-08 ***
## log10(length_av)  0.928475   0.091351  10.164 7.636e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.05883 on 14 degrees of freedom
## Multiple R-squared: 0.8807,  Adjusted R-squared: 0.8721 
## F-statistic: 103.3 on 1 and 14 DF,  p-value: 7.636e-08
#save coefficients of fits as object
cc_pgls_evo0 <- data.frame(coef(pgls_evo0))

Here is the model fit with lambda fixed at 1 (highest phylogenetic signal in model residuals). *Note that when we used a ML estimate of lambda, this was the value used, though the lambda estimate was not significant.

#fit model for eye diameter ~ body length (lambda = 1)
pgls_evo1 <- pgls(log10(eye_av) ~ log10(length_av), 
               data = species.comp, 
               lambda = 1, #set lambda to 1
               bounds=list(lambda=c(1,1)),
               param.CI = 0.95)

#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_evo1)

par(mfrow = c(1, 1))

#parameter estimates
summary(pgls_evo1)
## 
## Call:
## pgls(formula = log10(eye_av) ~ log10(length_av), data = species.comp, 
##     lambda = 1, param.CI = 0.95, bounds = list(lambda = c(1, 
##         1)))
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.100464 -0.041912  0.009141  0.044272  0.080878 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [Fix]  : 1.000
## delta  [Fix]  : 1.000
## 
## Coefficients:
##                   Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)      -1.288069   0.147956 -8.7057 5.051e-07 ***
## log10(length_av)  0.822852   0.095634  8.6042 5.810e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.05855 on 14 degrees of freedom
## Multiple R-squared: 0.841,   Adjusted R-squared: 0.8296 
## F-statistic: 74.03 on 1 and 14 DF,  p-value: 5.81e-07
#save coefficients of fits as object
cc_pgls_evo1 <- data.frame(coef(pgls_evo1))

We can plot both estimates of the relationship between eye size and body size onto our data to see the range of possible fits depending on the true value of lambda.

# Make plot 
plot_evo_bounds <- ggplot(species, aes(x = Body_length, y = Eye_diameter, text = genus_species)) + 
  geom_point(size = 2, alpha = 0.6) + 
  scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
  scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  geom_abline(data = cc_pgls_evo0, aes(intercept = cc_pgls_evo0[1,1], slope = cc_pgls_evo0[2,1]), linetype = "solid") +
  geom_abline(data = cc_pgls_evo1, aes(intercept = cc_pgls_evo1[1,1], slope = cc_pgls_evo1[2,1]), linetype = "dashed")
  

# Interactive plot
ggplotly(plot_evo_bounds)

Eye diameter vs. body length across 16 species of sergesteds. Solid line indicates the PGLS fit when lambda = 0 and dashed line the fit when lambda = 1.

We can also plot the static allometry for each species onto the evolutionary allometry (across species means). Most of the species have very similar static/intraspecific slopes to the evolutionary/interspecific pattern, but those few small-bodied species are the strange ones here too.

# Make plot 
plot_evo_bounds2 <- ggplot(species, aes(x = Body_length, y = Eye_diameter, text = genus_species)) + 
  geom_point(size = 2, alpha = 0.6) + 
  scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
  scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  geom_segment(aes(x = limits$min[which(limits$genus_species == "Allosergestes_pectinatus")], xend = limits$max[which(limits$genus_species == "Allosergestes_pectinatus")], y = 10^(log10(limits$min[which(limits$genus_species == "Allosergestes_pectinatus")])*cc_species["Allosergestes_pectinatus", "slope"] + cc_species["Allosergestes_pectinatus", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Allosergestes_pectinatus")])*cc_species["Allosergestes_pectinatus", "slope"] + cc_species["Allosergestes_pectinatus", "elevation"])), color = "grey85") + #Allosergestes pectinatus adult scaling
  geom_segment(aes(x = limits$min[which(limits$genus_species == "Allosergestes_sargassi")], xend = limits$max[which(limits$genus_species == "Allosergestes_sargassi")], y = 10^(log10(limits$min[which(limits$genus_species == "Allosergestes_sargassi")])*cc_species["Allosergestes_sargassi", "slope"] + cc_species["Allosergestes_sargassi", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Allosergestes_sargassi")])*cc_species["Allosergestes_sargassi", "slope"] + cc_species["Allosergestes_sargassi", "elevation"])), color = "grey85") + #Allosergestes_sargassi adult scaling
  geom_segment(aes(x = limits$min[which(limits$genus_species == "Challengerosergia_hansjacobi")], xend = limits$max[which(limits$genus_species == "Challengerosergia_hansjacobi")], y = 10^(log10(limits$min[which(limits$genus_species == "Challengerosergia_hansjacobi")])*cc_species["Challengerosergia_hansjacobi", "slope"] + cc_species["Challengerosergia_hansjacobi", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Challengerosergia_hansjacobi")])*cc_species["Challengerosergia_hansjacobi", "slope"] + cc_species["Challengerosergia_hansjacobi", "elevation"])), color = "grey85") + #Challengerosergia_hansjacobi adult scaling
  geom_segment(aes(x = limits$min[which(limits$genus_species == "Challengerosergia_talismani")], xend = limits$max[which(limits$genus_species == "Challengerosergia_talismani")], y = 10^(log10(limits$min[which(limits$genus_species == "Challengerosergia_talismani")])*cc_species["Challengerosergia_talismani", "slope"] + cc_species["Challengerosergia_talismani", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Challengerosergia_talismani")])*cc_species["Challengerosergia_talismani", "slope"] + cc_species["Challengerosergia_talismani", "elevation"])), color = "grey85") + #Challengerosergia_talismani adult scaling
  geom_segment(aes(x = limits$min[which(limits$genus_species == "Deosergestes_corniculum")], xend = limits$max[which(limits$genus_species == "Deosergestes_corniculum")], y = 10^(log10(limits$min[which(limits$genus_species == "Deosergestes_corniculum")])*cc_species["Deosergestes_corniculum", "slope"] + cc_species["Deosergestes_corniculum", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Deosergestes_corniculum")])*cc_species["Deosergestes_corniculum", "slope"] + cc_species["Deosergestes_corniculum", "elevation"])), color = "grey85") + #Deosergestes_corniculum adult scaling
  geom_segment(aes(x = limits$min[which(limits$genus_species == "Deosergestes_henseni")], xend = limits$max[which(limits$genus_species == "Deosergestes_henseni")], y = 10^(log10(limits$min[which(limits$genus_species == "Deosergestes_henseni")])*cc_species["Deosergestes_henseni", "slope"] + cc_species["Deosergestes_henseni", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Deosergestes_henseni")])*cc_species["Deosergestes_henseni", "slope"] + cc_species["Deosergestes_henseni", "elevation"])), color = "grey85") + #Deosergestes henseni adult scaling
  geom_segment(aes(x = limits$min[which(limits$genus_species == "Gardinerosergia_splendens")], xend = limits$max[which(limits$genus_species == "Gardinerosergia_splendens")], y = 10^(log10(limits$min[which(limits$genus_species == "Gardinerosergia_splendens")])*cc_species["Gardinerosergia_splendens", "slope"] + cc_species["Gardinerosergia_splendens", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Gardinerosergia_splendens")])*cc_species["Gardinerosergia_splendens", "slope"] + cc_species["Gardinerosergia_splendens", "elevation"])), color = "grey85") + #Gardinerosergia splendens adult scaling
  geom_segment(aes(x = limits$min[which(limits$genus_species == "Parasergestes_armatus")], xend = limits$max[which(limits$genus_species == "Parasergestes_armatus")], y = 10^(log10(limits$min[which(limits$genus_species == "Parasergestes_armatus")])*cc_species["Parasergestes_armatus", "slope"] + cc_species["Parasergestes_armatus", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Parasergestes_armatus")])*cc_species["Parasergestes_armatus", "slope"] + cc_species["Parasergestes_armatus", "elevation"])), color = "grey85") + #Parasergestes armatus adult scaling
  geom_segment(aes(x = limits$min[which(limits$genus_species == "Parasergestes_vigilax")], xend = limits$max[which(limits$genus_species == "Parasergestes_vigilax")], y = 10^(log10(limits$min[which(limits$genus_species == "Parasergestes_vigilax")])*cc_species["Parasergestes_vigilax", "slope"] + cc_species["Parasergestes_vigilax", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Parasergestes_vigilax")])*cc_species["Parasergestes_vigilax", "slope"] + cc_species["Parasergestes_vigilax", "elevation"])), color = "grey85") + #Parasergestes vigilax adult scaling
  geom_segment(aes(x = limits$min[which(limits$genus_species == "Phorcosergia_grandis")], xend = limits$max[which(limits$genus_species == "Phorcosergia_grandis")], y = 10^(log10(limits$min[which(limits$genus_species == "Phorcosergia_grandis")])*cc_species["Phorcosergia_grandis", "slope"] + cc_species["Phorcosergia_grandis", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Phorcosergia_grandis")])*cc_species["Phorcosergia_grandis", "slope"] + cc_species["Phorcosergia_grandis", "elevation"])), color = "grey85") + #Phorcosergia grandis adult scaling
  geom_segment(aes(x = limits$min[which(limits$genus_species == "Robustasergia_robusta")], xend = limits$max[which(limits$genus_species == "Robustasergia_robusta")], y = 10^(log10(limits$min[which(limits$genus_species == "Robustasergia_robusta")])*cc_species["Robustasergia_robusta", "slope"] + cc_species["Robustasergia_robusta", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Robustasergia_robusta")])*cc_species["Robustasergia_robusta", "slope"] + cc_species["Robustasergia_robusta", "elevation"])), color = "grey85") + #Robustasergia robusta adult scaling
  geom_segment(aes(x = limits$min[which(limits$genus_species == "Robustosergia_regalis")], xend = limits$max[which(limits$genus_species == "Robustosergia_regalis")], y = 10^(log10(limits$min[which(limits$genus_species == "Robustosergia_regalis")])*cc_species["Robustosergia_regalis", "slope"] + cc_species["Robustosergia_regalis", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Robustosergia_regalis")])*cc_species["Robustosergia_regalis", "slope"] + cc_species["Robustosergia_regalis", "elevation"])), color = "grey85") +
  geom_segment(aes(x = limits$min[which(limits$genus_species == "Sergestes_atlanticus")], xend = limits$max[which(limits$genus_species == "Sergestes_atlanticus")], y = 10^(log10(limits$min[which(limits$genus_species == "Sergestes_atlanticus")])*cc_species["Sergestes_atlanticus", "slope"] + cc_species["Sergestes_atlanticus", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Sergestes_atlanticus")])*cc_species["Sergestes_atlanticus", "slope"] + cc_species["Sergestes_atlanticus", "elevation"])), color = "grey85") +
  geom_segment(aes(x = limits$min[which(limits$genus_species == "Sergia_tenuiremis")], xend = limits$max[which(limits$genus_species == "Sergia_tenuiremis")], y = 10^(log10(limits$min[which(limits$genus_species == "Sergia_tenuiremis")])*cc_species["Sergia_tenuiremis", "slope"] + cc_species["Sergia_tenuiremis", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Sergia_tenuiremis")])*cc_species["Sergia_tenuiremis", "slope"] + cc_species["Sergia_tenuiremis", "elevation"])), color = "grey85") +
  geom_segment(aes(x = limits$min[which(limits$genus_species == "Eusergestes_arcticus")], xend = limits$max[which(limits$genus_species == "Eusergestes_arcticus")], y = 10^(log10(limits$min[which(limits$genus_species == "Eusergestes_arcticus")])*cc_species["Eusergestes_arcticus", "slope"] + cc_species["Eusergestes_arcticus", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Eusergestes_arcticus")])*cc_species["Eusergestes_arcticus", "slope"] + cc_species["Eusergestes_arcticus", "elevation"])), color = "grey85") +
  geom_abline(data = cc_pgls_evo0, aes(intercept = cc_pgls_evo0[1,1], slope = cc_pgls_evo0[2,1]), linetype = "solid") +
  geom_abline(data = cc_pgls_evo1, aes(intercept = cc_pgls_evo1[1,1], slope = cc_pgls_evo1[2,1]), linetype = "dashed")

# Interactive plot
plot_evo_bounds2

#export evolutionary allometry figures to pdf

#witout static allom
pdf("../Figures/evo_allom.pdf", width = 4, height = 2.5)
plot_evo_bounds
dev.off()
## quartz_off_screen 
##                 2
#with static allom
pdf("../Figures/evo_allom-static.pdf", width = 4, height = 2.5)
plot_evo_bounds2
dev.off()
## quartz_off_screen 
##                 2

Correlates of eye size in deep-sea sergestid shrimp

Next, we examine whether variation in eye size is correlated with ecological variables related to light regime and vision: depth and bioluminescent organs.

Phylogenetic distribution of eye size and ecology

Here, we examine how absolute and relative eye size vary across the sergestid phylogeny, and the corresponding distribution ecological variables. For visualizing relative eye size, we use the residuals of the PGLS fit for eye-body allometry, as these show eye size corrected for both body size and evolutionary relationships. (Here I have used the model where lambda = 1, as that was the ML estimate despite being non-significant.)

#Add residuals from PGLS fits to dataframe ----

#extract pgls residuals 
pglsres <- residuals(pgls_evo) 

#name residuals
colnames(pglsres) <- "pglsres" 

#add rownames to species dataframe
rownames(species) <- species$Binomial

#merge residuals with species data by rowname
species.phy <- merge(species, pglsres, by = "row.names") 

#make column converted to observed eye size as factor of PGLS fit
species.phy$eyefactor <- 10^species.phy$pglsres

# Prep phylogeny -----

# Make Genus_species tip labels to replace Binomial codes
sp.tips <- data.frame(old = species.phy$Binomial, new = species.phy$genus_species)

# replace tip labels in phylogeny
tree.plot <- sub.taxa.label(tree, sp.tips)

# Prep data ----

#change row names of species data
rownames(species.phy) <- species.phy$genus_species

#check that phylogeny and data match exactly
name.check(tree.plot, species.phy)

#resort trait dataset to the order of tree tip labels
species.phy <- species.phy[tree.plot$tip.label, ] 

#make trait vector for absolute eye size
aveye <- as.vector(species.phy$eye_av) 
names(aveye) <- species.phy$genus_species

#make trait vector for eye investment (PGLS residuals)
releye <- as.vector(species.phy$pglsres) #residuals of pgls
names(releye) <- species.phy$genus_species

# Phylogeny with eye size and investment -------

#color vector for photophores
col_ph <- c("lensed" = "palegreen",
            "none" = "honeydew",
            "pesta" = "cadetblue1",
            "unlensed" = "orchid1")

#plot tree with eye investment bars
plotTree.wBars(tree.plot, releye, 
               scale = 1.5, 
               tip.labels = TRUE,
               col = unname(col_ph[as.factor(species.phy$Organ)]),
               offset = 1,
               ftype = "bi",
               fsize = 0.7)

#add tip labels for absolute eye size
tiplabels(cex = 1.5*aveye,
          pch = 19, #shape of labels
          col = "mediumpurple",
          offset = 0)

#add legend for photophores
legend(x = "bottomleft", legend = c("lensed", "none", "pesta", "unlensed"), pch = 22, pt.cex= 2, pt.bg = col_ph, cex = 1, bty = "n", horiz = F)
Phylogeny of the 15 sertestids sampled depicting absolute eye diameter (purple circles scaled by eye diameter), relative eye investment (bars scaled by the residuals from the PGLS of eye diameter vs. body length, lambda = 1), and photophore complexity (color of bars).

Phylogeny of the 15 sertestids sampled depicting absolute eye diameter (purple circles scaled by eye diameter), relative eye investment (bars scaled by the residuals from the PGLS of eye diameter vs. body length, lambda = 1), and photophore complexity (color of bars).

#Export figure------
pdf("../Figures/phylo-photo1.pdf", width = 6, height = 4)

#plot tree with eye investment bars
plotTree.wBars(tree.plot, releye, 
               scale = 1.5, 
               tip.labels = TRUE,
               col = unname(col_ph[species.phy$Organ]),
               offset = 1,
               ftype = "bi",
               fsize = 0.7)

#add tip labels for absolute eye size
tiplabels(cex = 1.5*aveye,
          pch = 19, #shape of labels
          col = "mediumpurple",
          offset = 0)

#add legend for photophores
legend(x = "bottomleft", legend = c("lensed", "none", "pesta", "unlensed"), pch = 22, pt.cex= 2, pt.bg = col_ph, cex = 1, bty = "n", horiz = F)


dev.off()

We can also visualize day and night depth distributions.

# tree showing day and night mean depths across species

#fix node order upset by ladderization

#filter out internal nodes from second column of edge matrix
is_tip <- tree.plot$edge[,2] <= length(tree.plot$tip.label)
ordered_tips <- tree.plot$edge[is_tip, 2]

#extract tips in right order
#tree.plot$tip.label[ordered_tips]

#resort trait dataset to the order of tree tip labels
species.phy <- species.phy[tree.plot$tip.label[ordered_tips], ] 

#make trait vector for day depth av
day <- as.vector(species.phy$Day_Avg) 
names(day) <- species.phy$genus_species

#make trait vector for night depth av
night <- as.vector(species.phy$Night_Avg) 
names(night) <- species.phy$genus_species

#make trait vector for max reported depth
max <- as.vector(species.phy$Max_reported) 
names(max) <- species.phy$genus_species


#Plot tree with depth ranges
plot(tree.plot, 
     show.tip.label= TRUE, 
     font=3, 
     cex=0.9, 
     label.offset = 0.8, 
     edge.width=2)

#add tip labels for absolute eye size
tiplabels(cex = 1.5*aveye,
          pch = 19, #shape of labels
          col = "mediumpurple",
          offset = 0)

#add gray bar for difference in day/night depths (migration distance)
segments(1.05 + (day*0.0003), 1:length(day), 1.05 + (night*0.0003), 1:length(night), col="gray", lwd=8)

#add yellow point for day depth
segments(1.05 + (day*0.0003), 1:length(day), 1.05 + (day*0.0003), 1:length(day), col="orange", lwd=10)

#add blue point for night depth
segments(1.05 + (night*0.0003), 1:length(night), 1.05 + (night*0.0003), 1:length(night), col="royalblue2", lwd=10)

#add black point for max depth
segments(1.05 + (max*0.0003), 1:length(max), 1.05 + (max*0.0003), 1:length(max), col="black", lwd=7)

# export figure 
pdf("../Figures/phylo-depth.pdf", width = 8, height = 5.5)

#Plot tree with depth ranges
plot(tree.plot, 
     show.tip.label= TRUE, 
     font=3, 
     cex=0.9, 
     label.offset = 0.8, 
     edge.width=2)

#add tip labels for absolute eye size
tiplabels(cex = 1.5*aveye,
          pch = 19, #shape of labels
          col = "mediumpurple",
          offset = 0)

#add gray bar for difference in day/night depths (migration distance)
segments(1.05 + (day*0.0003), 1:length(day), 1.05 + (night*0.0003), 1:length(night), col="gray", lwd=8)

#add yellow point for day depth
segments(1.05 + (day*0.0003), 1:length(day), 1.05 + (day*0.0003), 1:length(day), col="orange", lwd=10)

#add blue point for night depth
segments(1.05 + (night*0.0003), 1:length(night), 1.05 + (night*0.0003), 1:length(night), col="royalblue2", lwd=10)

#add black point for max depth
segments(1.05 + (max*0.0003), 1:length(max), 1.05 + (max*0.0003), 1:length(max), col="black", lwd=7)

dev.off()

The above plot shows absolute eye size (purple dots scaled by eye diameter) and the average daytime (orange) and nighttime (blue) depths for each species in the phylogeny, as well as the maximum reported depth (black). Gray bars indicate vertical distance traveled during diel vertical migrations. Note that this plot is missing a scale bar, easier to add in Illustrator, but minimum depth at night (left/blue points) is 133m and maximum depth during the day (right/orange points) is 1750m. The maximum depth of occurrance (black points) across species is 2300m.

Photophore complexity

Does eye size correlate with photophore complexity?

Prediction: Eye size may be largest in animals with unlensed photophores, followed by lensed and organ of pesta types, as organs that are visually more difficult to resolve will require greater sensitivity to light (note: should this say longer focal lengths or something like that for higer acuity, not sensitivity?).

Preliminary plots (photophores)

Before incorporating phylogeny, we can inspect the distribution of our data with boxplots for absolute and relative eye size by photophore type.

#Boxplots for absolute eye size across photophores -----
plot_abs_photo <- ggplot(data = species.phy, aes(x = reorder(Organ, eye_av, FUN = mean), y = eye_av, text = genus_species)) + 
  geom_boxplot(notch = FALSE, outlier.shape = NA, alpha = 0.6) + #controls boxes
  stat_summary(fun.y = mean, colour = "black", geom = "point", shape = 18, size = 3, show_guide = FALSE) + #controls what stats shown
 geom_jitter(shape = 19, size = 1.3, alpha = 0.5, position = position_jitter(0.15)) +
  theme(text = element_text(size=14), panel.background = element_blank(), axis.line = element_line(colour = "black")) + #controls background
  xlab("Photophore complexity") +
  ylab("Eye diameter (mm)") +
  ggtitle("Absolute eye size")

ggplotly(plot_abs_photo)
#Boxplots for relative eye investment across photophores -----
plot_rel_photo <- ggplot(data = species.phy, aes(x = reorder(Organ, pglsres, FUN = mean), y = pglsres, text = genus_species)) + 
  geom_boxplot(notch = FALSE, outlier.shape = NA, alpha = 0.6) + #controls boxes
  stat_summary(fun.y = mean, colour = "black", geom = "point", shape = 18, size = 3, show_guide = FALSE) + #controls what stats shown
 geom_jitter(shape = 19, size = 1.3, alpha = 0.5, position = position_jitter(0.15)) +
  theme(text = element_text(size=14), panel.background = element_blank(), axis.line = element_line(colour = "black")) + #controls background
  xlab("Photophore complexity") +
  ylab("Eye ~ body residuals") +
  ggtitle("Relative eye investment")

ggplotly(plot_rel_photo) 

It looks like species with organs of pesta tend to have smaller eye diameters and smaller relative eye sizes than species with other types of photophores (note: we may want to combine the non-pesta organs together for power as sample size within groups is really small). We next test this in an evolutionary framework.

Test absolute eye size across photophores (PGLS)

This model is for eye diameter vs. photophore type. Lambda is estimated by maximum likelihood.

#run PGLS model using the ML estimate of lambda
pgls_abs_photo <- pgls(eye_av ~ Organ, 
               data = species.comp, 
               lambda = "ML",
               #uses Maximum Liklihood estimate of lambda
               param.CI = 0.95)

#evaluate model assumptions
par(mfrow = c(2,2)) #makes your plot window into 2x2 panels
plot(pgls_abs_photo) #plot the linear model

par(mfrow = c(1,1)) #set plotting window back to one plot

#main effects
anova(pgls_abs_photo)
## Analysis of Variance Table
## Sequential SS for pgls: lambda = 1.00, delta = 1.00, kappa = 1.00
## 
## Response: eye_av
##           Df  Sum Sq  Mean Sq F value  Pr(>F)  
## Organ      3 0.71605 0.238683  2.8559 0.08159 .
## Residuals 12 1.00292 0.083576                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#coefficients
summary(pgls_abs_photo)
## 
## Call:
## pgls(formula = eye_av ~ Organ, data = species.comp, lambda = "ML", 
##     param.CI = 0.95)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.55499 -0.18875 -0.03819  0.05956  0.27422 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 1.000
##    lower bound : 0.000, p = 0.21126
##    upper bound : 1.000, p = 1    
##    95.0% CI   : (NA, NA)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##               Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)    0.87918    0.12205  7.2035 1.081e-05 ***
## Organlensed    0.21111    0.26870  0.7857   0.44729    
## Organunlensed  0.56885    0.20812  2.7332   0.01816 *  
## Organnone      0.32053    0.28563  1.1222   0.28373    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2891 on 12 degrees of freedom
## Multiple R-squared: 0.4166,  Adjusted R-squared: 0.2707 
## F-statistic: 2.856 on 3 and 12 DF,  p-value: 0.08159
#Likelihood plot for Pagel's lambda 
plot(pgls.profile(pgls_abs_photo, "lambda"), 
     main = "Pagel's lambda", 
     sub = "Solid red line indicates estimate for lambda and broken red lines indcaite the 95% confidence interval")

This test shows no significant effect of photophore type (pesta, lensed, unlensed, none) on mean absolute eye size across species. The ML estimate of lambda in this model is 1, but this is not significant. So, next we fix lambda at its bounds (0 and 1) to see the effects on our results.

Test when lambda = 0

Absolute eye size vs. photophore type

#fit model for eye diameter ~ body length (lambda = 0)
pgls_abs_photo0 <- pgls(eye_av ~ Organ, 
               data = species.comp, 
               lambda = 0.00001, #set lambda to 0
               bounds=list(lambda=c(0.00001,0.00001)),
               param.CI = 0.95)

#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_abs_photo0)

par(mfrow = c(1, 1))

#main effects
anova(pgls_abs_photo0)
## Analysis of Variance Table
## Sequential SS for pgls: lambda = 0.00, delta = 1.00, kappa = 1.00
## 
## Response: eye_av
##           Df  Sum Sq Mean Sq F value  Pr(>F)  
## Organ      3 1.24698 0.41566  5.4302 0.01361 *
## Residuals 12 0.91855 0.07655                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#parameter estimates
summary(pgls_abs_photo0)
## 
## Call:
## pgls(formula = eye_av ~ Organ, data = species.comp, lambda = 1e-05, 
##     param.CI = 0.95, bounds = list(lambda = c(1e-05, 1e-05)))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.49463 -0.23850 -0.09265  0.12102  0.48200 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [Fix]  : 0.000
## delta  [Fix]  : 1.000
## 
## Coefficients:
##               Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)   0.829177   0.092223  8.9910 1.116e-06 ***
## Organlensed   0.319005   0.216283  1.4749  0.165975    
## Organunlensed 0.653861   0.166258  3.9328  0.001989 ** 
## Organnone     0.420822   0.291635  1.4430  0.174619    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2767 on 12 degrees of freedom
## Multiple R-squared: 0.5758,  Adjusted R-squared: 0.4698 
## F-statistic:  5.43 on 3 and 12 DF,  p-value: 0.01361

When lambda is set to 0, there is a significant effect of organ type on absolute eye size.

Test when lambda = 1

Absolute eye size vs. photophore type.

#fit model
pgls_abs_photo1 <- pgls(eye_av ~ Organ, 
               data = species.comp, 
               lambda = 1, #set lambda to 1
               bounds=list(lambda=c(1,1)),
               param.CI = 0.95)

#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_abs_photo1)

par(mfrow = c(1, 1))

#main effects
anova(pgls_abs_photo1)
## Analysis of Variance Table
## Sequential SS for pgls: lambda = 1.00, delta = 1.00, kappa = 1.00
## 
## Response: eye_av
##           Df  Sum Sq  Mean Sq F value  Pr(>F)  
## Organ      3 0.71605 0.238683  2.8559 0.08159 .
## Residuals 12 1.00292 0.083576                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#parameter estimates
summary(pgls_abs_photo1)
## 
## Call:
## pgls(formula = eye_av ~ Organ, data = species.comp, lambda = 1, 
##     param.CI = 0.95, bounds = list(lambda = c(1, 1)))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.55499 -0.18875 -0.03819  0.05956  0.27422 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [Fix]  : 1.000
## delta  [Fix]  : 1.000
## 
## Coefficients:
##               Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)    0.87918    0.12205  7.2035 1.081e-05 ***
## Organlensed    0.21111    0.26870  0.7857   0.44729    
## Organunlensed  0.56885    0.20812  2.7332   0.01816 *  
## Organnone      0.32053    0.28563  1.1222   0.28373    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2891 on 12 degrees of freedom
## Multiple R-squared: 0.4166,  Adjusted R-squared: 0.2707 
## F-statistic: 2.856 on 3 and 12 DF,  p-value: 0.08159

When lambda is set to 1, there is no significant effect of photophore type (pesta, lensed, unlensed, none) on mean absolute eye diameter across species.

Note here that though our overall test for significant effect of photophore type (the anova) is not significant (p = 0.08), the contrasts show one significant difference between the unlensed and the pesta states. The contrasts here are more post-hoc, so we shouldn’t really consider them given that the full model was not significant. Full model wins out. For the paper, I’d only report these contrasts for models with a significant effect of the ecological trait being tested.

Test relative eye size across photophores (PGLS)

Below, I run a PGLS model of log10(eye diameter) vs. log10(body length) with photophore complexity as a covariate.

#run PGLS model using the Maximum Liklihood estimate of lambda
pgls_photo <- pgls(log10(eye_av) ~ log10(length_av) + Organ, 
               data = species.comp, 
               lambda = "ML",
               #uses Maximum Liklihood estimate of lambda
               param.CI = 0.95)

#evaluate model assumptions
par(mfrow = c(2,2)) #makes your plot window into 2x2 panels
plot(pgls_photo) #plot the linear model

par(mfrow = c(1,1)) #set plotting window back to one plot

The model checks look pretty good here.

Next we can take a look at the overall significance of the main effects and sequential sums of squares. Here we can see that body length has a significant effect on eye size (expected, we already knew that), but we also find that photophore complexity/type has a significant effect on eye size when the effect of body size is accounted for (p = 0.01).

#anova
anova(pgls_photo)
## Analysis of Variance Table
## Sequential SS for pgls: lambda = 0.00, delta = 1.00, kappa = 1.00
## 
## Response: log10(eye_av)
##                  Df  Sum Sq Mean Sq  F value    Pr(>F)    
## log10(length_av)  1 0.35757 0.35757 213.2573 1.511e-08 ***
## Organ             3 0.03002 0.01001   5.9672   0.01144 *  
## Residuals        11 0.01844 0.00168                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Finally, we look at the PGLS coefficients as well as significance with respect to contrasts. Note that for a categorical variable, R will treat whichever factor comes first as the “control” and will test other factors for significant differences relative to that “control”. Here, the comparison state (“control”) for photophore complexity is set to “pesta” (estimate listed under (Intercept) - that’s why you don’t see that state as a comparison line) and the other states show the difference in the estimated intercept for that state compared to the pesta state, and whether that difference is significant. You can see that “unlensed” and “lensed” both show significant differences from “pesta”, and the intercept for both is higher, indicating larger mean relative eye sizes in those groups.

#summary
summary(pgls_photo)
## 
## Call:
## pgls(formula = log10(eye_av) ~ log10(length_av) + Organ, data = species.comp, 
##     lambda = "ML", param.CI = 0.95)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.083912 -0.026000 -0.007001  0.014353  0.041549 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 0.000
##    lower bound : 0.000, p = 1    
##    upper bound : 1.000, p = 0.099077
##    95.0% CI   : (NA, NA)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##                   Estimate Std. Error  t value  Pr(>|t|)    
## (Intercept)      -1.229881   0.114139 -10.7753 3.489e-07 ***
## log10(length_av)  0.769679   0.077417   9.9420 7.835e-07 ***
## Organlensed       0.099594   0.032644   3.0509  0.011029 *  
## Organunlensed     0.108865   0.029393   3.7038  0.003478 ** 
## Organnone         0.051706   0.045672   1.1321  0.281667    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.04095 on 11 degrees of freedom
## Multiple R-squared: 0.9546,  Adjusted R-squared: 0.9381 
## F-statistic: 57.79 on 4 and 11 DF,  p-value: 2.576e-07

The output also includes the maximum-likelihood estimation of lambda, which represents the phylogenetic signal for PGLS model residuals. Here, the estimate for lambda is 0.00, which means that there is no phylogenetic signal. This doesn’t mean we shouldn’t include phylogeny in the model (no reason not to), but it does mean that it’s the same result we would get if we assessed this without including phylogeny. This also doesn’t mean that there is not phylogenetic signal in the individual traits that we put into the model; this is only lambda for the model residuals, not the traits themselves. However, as mentioned above, when you have <20 species the power to estimate lambda is not there, so we may want to consider running these models twice, one forcing lambda to 1 and another forcing lambda to 0, so that we can see how that affects our findings. That could help us in our interpretation (because photophore type is tightly matched to clades, my expectation is that if we force lambda to 0.5 or 1, the relationship may not remain significant.)

#Likelihood plot for Pagel's lambda 
plot(pgls.profile(pgls_photo, "lambda"), 
     main = "Pagel's lambda", 
     sub = "Solid red line indicates estimate for lambda and broken red lines indcaite the 95% confidence interval")

With sample sizes too small to generate a robust estimate of lambda, one option is to run models with the lowest and possible bounds of lambda (0 and 1), so that we can see what the range of possible parameter estimates are depending on the phylogenetic signal of the model residuals. Below, we fit these two models so that we can compare them.

Test when lambda = 0

Here is the model fit with lambda fixed at 0 (no phylogenetic signal in model residuals). *Note that when we used a ML estimate of lambda, this was the value used, though the lambda estimate was not significant.

#fit model for eye diameter ~ body length (lambda = 0)
pgls_photo0 <- pgls(log10(eye_av) ~ log10(length_av) + Organ, 
               data = species.comp, 
               lambda = 0.00001, #set lambda to 0
               bounds=list(lambda=c(0.00001,0.00001)),
               param.CI = 0.95)

#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_photo0)

par(mfrow = c(1, 1))

#parameter estimates
summary(pgls_photo0)
## 
## Call:
## pgls(formula = log10(eye_av) ~ log10(length_av) + Organ, data = species.comp, 
##     lambda = 1e-05, param.CI = 0.95, bounds = list(lambda = c(1e-05, 
##         1e-05)))
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.083912 -0.026001 -0.007001  0.014353  0.041549 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [Fix]  : 0.000
## delta  [Fix]  : 1.000
## 
## Coefficients:
##                   Estimate Std. Error  t value  Pr(>|t|)    
## (Intercept)      -1.229880   0.114139 -10.7753 3.489e-07 ***
## log10(length_av)  0.769678   0.077417   9.9420 7.835e-07 ***
## Organlensed       0.099594   0.032644   3.0509  0.011029 *  
## Organunlensed     0.108865   0.029393   3.7038  0.003478 ** 
## Organnone         0.051706   0.045673   1.1321  0.281667    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.04095 on 11 degrees of freedom
## Multiple R-squared: 0.9546,  Adjusted R-squared: 0.9381 
## F-statistic: 57.79 on 4 and 11 DF,  p-value: 2.576e-07
anova(pgls_photo0)
## Analysis of Variance Table
## Sequential SS for pgls: lambda = 0.00, delta = 1.00, kappa = 1.00
## 
## Response: log10(eye_av)
##                  Df  Sum Sq Mean Sq  F value    Pr(>F)    
## log10(length_av)  1 0.35757 0.35757 213.2555 1.511e-08 ***
## Organ             3 0.03002 0.01001   5.9671   0.01144 *  
## Residuals        11 0.01844 0.00168                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#save coefficients of fits as object
cc_pgls_photo0 <- data.frame(coef(pgls_photo0))

Here is the model fit with lambda fixed at 1 (highest phylogenetic signal in model residuals).

#fit model for eye diameter ~ body length (lambda = 1)
pgls_photo1 <- pgls(log10(eye_av) ~ log10(length_av) + Organ, 
               data = species.comp, 
               lambda = 1, #set lambda to 1
               bounds=list(lambda=c(1,1)),
               param.CI = 0.95)

#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_photo1)

par(mfrow = c(1, 1))

#model outputs
anova(pgls_photo1)
## Analysis of Variance Table
## Sequential SS for pgls: lambda = 1.00, delta = 1.00, kappa = 1.00
## 
## Response: log10(eye_av)
##                  Df   Sum Sq  Mean Sq  F value    Pr(>F)    
## log10(length_av)  1 0.253779 0.253779 106.0654 5.504e-07 ***
## Organ             3 0.021672 0.007224   3.0193   0.07574 .  
## Residuals        11 0.026319 0.002393                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(pgls_photo1)
## 
## Call:
## pgls(formula = log10(eye_av) ~ log10(length_av) + Organ, data = species.comp, 
##     lambda = 1, param.CI = 0.95, bounds = list(lambda = c(1, 
##         1)))
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.067903 -0.023288 -0.004012  0.039617  0.064583 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [Fix]  : 1.000
## delta  [Fix]  : 1.000
## 
## Coefficients:
##                   Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)      -1.185672   0.135069 -8.7783 2.674e-06 ***
## log10(length_av)  0.740092   0.089218  8.2953 4.620e-06 ***
## Organlensed       0.102270   0.045518  2.2468   0.04615 *  
## Organunlensed     0.107974   0.038334  2.8167   0.01677 *  
## Organnone         0.066800   0.049589  1.3471   0.20505    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.04891 on 11 degrees of freedom
## Multiple R-squared: 0.9128,  Adjusted R-squared: 0.8811 
## F-statistic: 28.78 on 4 and 11 DF,  p-value: 8.972e-06

Photophores have a significant effect on eye size when lambda is set to 0 (no phylogenetic signal in model residuals). However, we lose that significance when lambda is set to 1 (high phylogenetic signal in model residuals). While not significant here, it is close (p = 0.08). We can’t be totally confident in our finding, as it depends on the true value of lambda (which we can’t estimate with our sample size of species). However, since all patterns we see seem to be driven by organ of pesta, maybe we could try combining categories into organ of pesta or not, for reduced states and more power? I have done that next.

Photophore complexity (pesta vs. not - reduced states)

Here, I’ve reduced the number of states to the following: bioluminescent with organ of pesta present (“pesta”), bioluminescent with organ of pesta not present (“nonpesta”), bioluminescence not present (“none”). Then I’ve re-run the models with lambda set to 0 and 1 to see the binned range of potential findings depending on the true value of lambda.

#recategorize photophore bins
species2 <- species %>%
  mutate(photo2 = recode(Organ, "lensed" = "nonpesta", "unlensed" = "nonpesta")) %>%
  mutate(photo2 = factor(photo2, levels = c("pesta","nonpesta","none")))

#use caper function to combine phylogeny and data into one object (this function also matches species names in tree and dataset)
species.comp2 <- comparative.data(data = species2, phy = tree, names.col = "Binomial", vcv = TRUE, na.omit = FALSE, warn.dropped = TRUE)

#check for dropped tips or dropped species
species.comp2$dropped$tips #phylogeny
species.comp2$dropped$unmatched.rows #dataset

Data visualization

First we can look at the distribution of data again in these new states with a phylogeny visualization and with boxplots.

# Phylogeny with eye size and investment -------

#color palette for photophore groups
col_photo2 <- c("pesta" = "#7b3294",
                "lensed" = "#008837",
                "unlensed" = "#a6dba0",
                "none" = "#A6ACAF")

#plot tree with eye investment bars
plotTree.wBars(tree.plot, releye, 
               scale = 1.5, 
               tip.labels = TRUE,
               col = "gray",
               offset = 1,
               ftype = "bi",
               fsize = 0.7)

#add tip labels for absolute eye size
tiplabels(cex = 1.5*aveye,
          pch = 19, #shape of labels
          col = unname(col_photo2[species2$Organ]),
          offset = 0)

#add legend for photophores
legend(x = "bottomleft", legend = c("pesta", "lensed", "unlensed", "none"), pch = 22, pt.cex= 2, pt.bg = col_photo2, cex = 1, bty = "n", horiz = F)
Phylogeny of the 15 sertestids sampled depicting absolute eye diameter (circles, colored by photophore state) and relative eye investment (bars scaled by the residuals from the PGLS of eye diameter vs. body length, lambda = 1).

Phylogeny of the 15 sertestids sampled depicting absolute eye diameter (circles, colored by photophore state) and relative eye investment (bars scaled by the residuals from the PGLS of eye diameter vs. body length, lambda = 1).

#Export figure------
pdf("../Figures/phylo-photo2.pdf", width = 6, height = 4)

#plot tree with eye investment bars
plotTree.wBars(tree.plot, releye, 
               scale = 1.5, 
               tip.labels = TRUE,
               col = "gray",
               offset = 1,
               ftype = "bi",
               fsize = 0.7)

#add tip labels for absolute eye size
tiplabels(cex = 1.5*aveye,
          pch = 19, #shape of labels
          col = unname(col_photo2[species2$Organ]),
          offset = 0)

#add legend for photophores
legend(x = "bottomleft", legend = c("pesta", "lensed", "unlensed", "none"), pch = 22, pt.cex= 2, pt.bg = col_photo2, cex = 1, bty = "n", horiz = F)

dev.off()
sh_photo2 <- c("lensed" = 17, 
                "none" = 18,
                "pesta" = 19,
                "unlensed" = 15)

#Boxplots for absolute eye size across photophores -----
plot_abs_photo2 <- ggplot(data = species2, aes(x = factor(photo2, levels=c("pesta", "nonpesta","none")), y = eye_av)) + 
  geom_boxplot(notch = FALSE, outlier.shape = NA, alpha = 0.6) + #controls boxes
  stat_summary(fun.y = mean, colour = "black", geom = "point", shape = 18, size = 3, show_guide = FALSE) + #controls what stats shown
 geom_jitter(aes(x = photo2, y = eye_av, color = Organ, shape = Organ), size = 3, alpha = 0.9, position = position_jitter(0.15)) +
  theme(text = element_text(size=14), panel.background = element_blank(), axis.line = element_line(colour = "black"), legend.key = element_blank()) + #controls background
  scale_color_manual(values = col_photo2, name = "Photophore type", 
                     breaks = c("pesta", "lensed", "unlensed", "none")) +
  scale_shape_manual(values = sh_photo2, name = "Photophore type", 
                     breaks = c("pesta", "lensed", "unlensed", "none")) +
  xlab("Photophore complexity") +
  ylab("Eye diameter (mm)") +
  ggtitle("Absolute eye size")

plot_abs_photo2

#Boxplots for relative eye size across photophores -----
plot_rel_photo2 <- ggplot(data = species2, aes(x = factor(photo2, levels=c("pesta", "nonpesta","none")), y = eye_av/length_av)) + 
  geom_boxplot(notch = FALSE, outlier.shape = NA, alpha = 0.6) + #controls boxes
  stat_summary(fun.y = mean, colour = "black", geom = "point", shape = 18, size = 3, show_guide = FALSE) + #controls what stats shown
 geom_jitter(aes(x = photo2, y = eye_av/length_av, color = Organ, shape = Organ), size = 3, alpha = 0.9, position = position_jitter(0.15)) +
  theme(text = element_text(size=14), panel.background = element_blank(), axis.line = element_line(colour = "black"), legend.key = element_blank()) + #controls background
  scale_color_manual(values = col_photo2, name = "Photophore type", 
                     breaks = c("pesta", "lensed", "unlensed", "none")) +
  scale_shape_manual(values = sh_photo2, name = "Photophore type", 
                     breaks = c("pesta", "lensed", "unlensed", "none")) +
  xlab("Photophore complexity") +
  ylab("Relative eye size (eye/body)") +
  ggtitle("Relative eye size")

plot_rel_photo2

#Boxplots for residual eye investment across photophores -----
plot_inv_photo2 <- ggplot(data = species2, aes(x = factor(photo2, levels=c("pesta", "nonpesta","none")), y = pglsres)) + 
  geom_boxplot(notch = FALSE, outlier.shape = NA, alpha = 0.6) + #controls boxes
  stat_summary(fun.y = mean, colour = "black", geom = "point", shape = 18, size = 3, show_guide = FALSE) + #controls what stats shown
 geom_jitter(aes(x = photo2, y = pglsres, color = Organ, shape = Organ), size = 3, alpha = 0.9, position = position_jitter(0.15)) +
  theme(text = element_text(size=14), panel.background = element_blank(), axis.line = element_line(colour = "black"), legend.key = element_blank()) + #controls background
  scale_color_manual(values = col_photo2, name = "Photophore type", 
                     breaks = c("pesta", "lensed", "unlensed", "none")) +
  scale_shape_manual(values = sh_photo2, name = "Photophore type", 
                     breaks = c("pesta", "lensed", "unlensed", "none")) +
  xlab("Photophore complexity") +
  ylab("Eye ~ body residuals") +
  ggtitle("Residual eye investment")

plot_inv_photo2

# Make boxplot figure to export ----

#name panels
fig.a <- plot_abs_photo2 + theme(legend.position = "none") + ggtitle("")  + xlab("")
fig.b <- plot_rel_photo2 + theme(legend.position = "none") + ggtitle("") + xlab("")
fig.c <- plot_inv_photo2 + theme(legend.position = "none") + ggtitle("") + xlab("")

#arrange plots in panels
plots <- plot_grid(fig.a, fig.b, fig.c, #list of plots to arrange in grid
           align = 'vh', #align horizontally and vertically
           axis = 'lb',
           labels = c("A", "B", "C"), #panel labels for figure
           hjust = -1, #adjustment for panel labels
           nrow = 1) #number of rows in grids

# extract legend from Rana temporaria figure
leg <- get_legend(fig.a + theme(legend.position="right"))


#export figure
pdf("../Figures/boxplots-photo.pdf", width = 10, height = 4)
plot_grid(plots, leg, nrow = 1, rel_widths = c(1, .2))
dev.off()
## quartz_off_screen 
##                 2

Test absolute eye size vs. organ type (pesta, nonpesta, none)

This is a model of absolute eye diameter vs. organ type (PGLS). Again, we will first run a model that estimates lambda by ML.

#fit model 
pgls_photo_bin_abs <- pgls(eye_av ~ photo2, 
               data = species.comp2, 
               lambda = "ML",
               param.CI = 0.95)

#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_photo_bin_abs)

par(mfrow = c(1, 1))

#main effects
anova(pgls_photo_bin_abs)
## Analysis of Variance Table
## Sequential SS for pgls: lambda = 1.00, delta = 1.00, kappa = 1.00
## 
## Response: eye_av
##           Df  Sum Sq  Mean Sq F value Pr(>F)
## photo2     2 0.50453 0.252266  2.7004 0.1045
## Residuals 13 1.21443 0.093418
#parameter estimates
summary(pgls_photo_bin_abs)
## 
## Call:
## pgls(formula = eye_av ~ photo2, data = species.comp2, lambda = "ML", 
##     param.CI = 0.95)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.56167 -0.26508 -0.06907  0.06838  0.29391 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 1.000
##    lower bound : 0.000, p = 0.33913
##    upper bound : 1.000, p = 1    
##    95.0% CI   : (NA, NA)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##                Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)     0.87918    0.12904  6.8135 1.237e-05 ***
## photo2nonpesta  0.49215    0.21405  2.2992   0.03871 *  
## photo2none      0.33982    0.30171  1.1263   0.28039    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3056 on 13 degrees of freedom
## Multiple R-squared: 0.2935,  Adjusted R-squared: 0.1848 
## F-statistic:   2.7 on 2 and 13 DF,  p-value: 0.1045

When lambda is fitted by ML, the estimate for lambda is 1; however, this lambda estimate is not significant. Again, we can run two models bounding the extreme possibilities of lambda to see the range of possible results.

Test with lambda = 0

Absolute eye size vs. organ type (pesta, non-pesta, none)

#fit model for eye diameter ~ body length (lambda = 0)
pgls_photo0_bin_abs <- pgls(eye_av ~ photo2, 
               data = species.comp2, 
               lambda = 0.00001, #set lambda to 0
               bounds=list(lambda=c(0.00001,0.00001)),
               param.CI = 0.95)

#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_photo0_bin_abs)

par(mfrow = c(1, 1))

#main effects
anova(pgls_photo0_bin_abs)
## Analysis of Variance Table
## Sequential SS for pgls: lambda = 0.00, delta = 1.00, kappa = 1.00
## 
## Response: eye_av
##           Df Sum Sq Mean Sq F value  Pr(>F)  
## photo2     2 1.0975 0.54874  6.6791 0.01011 *
## Residuals 13 1.0680 0.08216                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#parameter estimates
summary(pgls_photo0_bin_abs)
## 
## Call:
## pgls(formula = eye_av ~ photo2, data = species.comp2, lambda = 1e-05, 
##     param.CI = 0.95, bounds = list(lambda = c(1e-05, 1e-05)))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.57427 -0.26319 -0.03726  0.13291  0.39016 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [Fix]  : 0.000
## delta  [Fix]  : 1.000
## 
## Coefficients:
##                Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)    0.829177   0.095544  8.6785 9.072e-07 ***
## photo2nonpesta 0.542243   0.151069  3.5894  0.003298 ** 
## photo2none     0.420823   0.302137  1.3928  0.187036    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2866 on 13 degrees of freedom
## Multiple R-squared: 0.5068,  Adjusted R-squared: 0.4309 
## F-statistic: 6.679 on 2 and 13 DF,  p-value: 0.01011

When lambda = 0, there is a significant effect of photophore type on absolute eye size, and the contrasts show that species with organs of pesta have significantly smaller mean relative eye sizes than those with other types of organ.

Test with lambda = 1

#fit model for eye diameter ~ body length (lambda = 1)
pgls_photo1_bin_abs <- pgls(eye_av ~ photo2, 
               data = species.comp2, 
               lambda = 1, #set lambda to 1
               bounds=list(lambda=c(1,1)),
               param.CI = 0.95)

#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_photo1_bin_abs)

par(mfrow = c(1, 1))

#main effects
anova(pgls_photo1_bin_abs)
## Analysis of Variance Table
## Sequential SS for pgls: lambda = 1.00, delta = 1.00, kappa = 1.00
## 
## Response: eye_av
##           Df  Sum Sq  Mean Sq F value Pr(>F)
## photo2     2 0.50453 0.252266  2.7004 0.1045
## Residuals 13 1.21443 0.093418
#parameter estimates
summary(pgls_photo1_bin_abs)
## 
## Call:
## pgls(formula = eye_av ~ photo2, data = species.comp2, lambda = 1, 
##     param.CI = 0.95, bounds = list(lambda = c(1, 1)))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.56167 -0.26508 -0.06907  0.06838  0.29391 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [Fix]  : 1.000
## delta  [Fix]  : 1.000
## 
## Coefficients:
##                Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)     0.87918    0.12904  6.8135 1.237e-05 ***
## photo2nonpesta  0.49215    0.21405  2.2992   0.03871 *  
## photo2none      0.33982    0.30171  1.1263   0.28039    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3056 on 13 degrees of freedom
## Multiple R-squared: 0.2935,  Adjusted R-squared: 0.1848 
## F-statistic:   2.7 on 2 and 13 DF,  p-value: 0.1045

When lambda = 1, there is no significant effect of photophore type on absolute eye size (this kind of makes sense, as all the species with the organs of pesta are one clade with one origin).

Test relative eye size vs. organ type

Here is a model of log10(eye diameter) vs. log10(body length) + photophore type (pesta, nonpesta, none), with lambda fitted by ML.

#fit model 
pgls_photo_bin <- pgls(log10(eye_av) ~ log10(length_av) + photo2, 
               data = species.comp2, 
               lambda = "ML",
               param.CI = 0.95)

#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_photo_bin)

par(mfrow = c(1, 1))

#main effects
anova(pgls_photo_bin)
## Analysis of Variance Table
## Sequential SS for pgls: lambda = 0.00, delta = 1.00, kappa = 1.00
## 
## Response: log10(eye_av)
##                  Df  Sum Sq Mean Sq  F value    Pr(>F)    
## log10(length_av)  1 0.35757 0.35757 231.3066 3.318e-09 ***
## photo2            2 0.02991 0.01495   9.6738  0.003147 ** 
## Residuals        12 0.01855 0.00155                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#parameter estimates
summary(pgls_photo_bin)
## 
## Call:
## pgls(formula = log10(eye_av) ~ log10(length_av) + photo2, data = species.comp2, 
##     lambda = "ML", param.CI = 0.95)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.084884 -0.027957 -0.003879  0.013107  0.042196 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 0.000
##    lower bound : 0.000, p = 1    
##    upper bound : 1.000, p = 0.10307
##    95.0% CI   : (NA, NA)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##                   Estimate Std. Error  t value  Pr(>|t|)    
## (Intercept)      -1.237405   0.105786 -11.6973 6.425e-08 ***
## log10(length_av)  0.774819   0.071712  10.8045 1.543e-07 ***
## photo2nonpesta    0.104921   0.023898   4.3904 0.0008801 ***
## photo2none        0.050714   0.043692   1.1607 0.2683204    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.03932 on 12 degrees of freedom
## Multiple R-squared: 0.9543,  Adjusted R-squared: 0.9429 
## F-statistic: 83.55 on 3 and 12 DF,  p-value: 2.614e-08

When lambda is fitted by ML, the estimate for lambda is 0, though this estimate is not significant. So, we fit at both extremes of lambda to see the range of possible results.

Test where lambda = 0

Here is the model fit with lambda fixed at 0 (no phylogenetic signal in model residuals).

#fit model for eye diameter ~ body length (lambda = 0)
pgls_photo0_bin <- pgls(log10(eye_av) ~ log10(length_av) + photo2, 
               data = species.comp2, 
               lambda = 0.00001, #set lambda to 0
               bounds=list(lambda=c(0.00001,0.00001)),
               param.CI = 0.95)

#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_photo0_bin)

par(mfrow = c(1, 1))

#main effects
anova(pgls_photo0_bin)
## Analysis of Variance Table
## Sequential SS for pgls: lambda = 0.00, delta = 1.00, kappa = 1.00
## 
## Response: log10(eye_av)
##                  Df  Sum Sq Mean Sq  F value    Pr(>F)    
## log10(length_av)  1 0.35757 0.35757 231.3046 3.318e-09 ***
## photo2            2 0.02991 0.01495   9.6737  0.003147 ** 
## Residuals        12 0.01855 0.00155                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#parameter estimates
summary(pgls_photo0_bin)
## 
## Call:
## pgls(formula = log10(eye_av) ~ log10(length_av) + photo2, data = species.comp2, 
##     lambda = 1e-05, param.CI = 0.95, bounds = list(lambda = c(1e-05, 
##         1e-05)))
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.084884 -0.027957 -0.003879  0.013107  0.042196 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [Fix]  : 0.000
## delta  [Fix]  : 1.000
## 
## Coefficients:
##                   Estimate Std. Error  t value  Pr(>|t|)    
## (Intercept)      -1.237405   0.105786 -11.6973 6.425e-08 ***
## log10(length_av)  0.774819   0.071712  10.8045 1.543e-07 ***
## photo2nonpesta    0.104921   0.023898   4.3904 0.0008801 ***
## photo2none        0.050714   0.043692   1.1607 0.2683199    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.03932 on 12 degrees of freedom
## Multiple R-squared: 0.9543,  Adjusted R-squared: 0.9429 
## F-statistic: 83.55 on 3 and 12 DF,  p-value: 2.614e-08

Photophore type has a significant effect on relative eye size when lambda = 0.

Here is the model fit with lambda fixed at 1 (highest phylogenetic signal in model residuals).

#fit model for eye diameter ~ body length (lambda = 1)
pgls_photo1_bin <- pgls(log10(eye_av) ~ log10(length_av) + photo2, 
               data = species.comp2, 
               lambda = 1, #set lambda to 1
               bounds=list(lambda=c(1,1)),
               param.CI = 0.95)

#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_photo1_bin)

par(mfrow = c(1, 1))

#main effects
anova(pgls_photo1_bin)
## Analysis of Variance Table
## Sequential SS for pgls: lambda = 1.00, delta = 1.00, kappa = 1.00
## 
## Response: log10(eye_av)
##                  Df   Sum Sq  Mean Sq  F value    Pr(>F)    
## log10(length_av)  1 0.253779 0.253779 115.4961 1.636e-07 ***
## photo2            2 0.021624 0.010812   4.9206   0.02751 *  
## Residuals        12 0.026368 0.002197                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#coefficients and contrasts
summary(pgls_photo1_bin)
## 
## Call:
## pgls(formula = log10(eye_av) ~ log10(length_av) + photo2, data = species.comp2, 
##     lambda = 1, param.CI = 0.95, bounds = list(lambda = c(1, 
##         1)))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.06902 -0.02289 -0.00536  0.03717  0.06538 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [Fix]  : 1.000
## delta  [Fix]  : 1.000
## 
## Coefficients:
##                   Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)      -1.191766   0.122725 -9.7109 4.912e-07 ***
## log10(length_av)  0.744165   0.080955  9.1923 8.827e-07 ***
## photo2nonpesta    0.106186   0.034696  3.0605   0.00989 ** 
## photo2none        0.066568   0.047496  1.4016   0.18638    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.04688 on 12 degrees of freedom
## Multiple R-squared: 0.9126,  Adjusted R-squared: 0.8908 
## F-statistic: 41.78 on 3 and 12 DF,  p-value: 1.255e-06

Here when we bin into 3 groups for photophores (organ of pesta, no organ of pesta, no photophores at all), we find that organ type has a significant effect on eye size in both models (lambda = 1 and lambda = 0). This is a robust result, as whatever the true value of lambda is, the association should hold. The only caveat I will add to this is that since we only have a single evolutionary transition from organ of pesta to non organ of pesta represented in our dataset, the association would be more convincing if found later across larger samples/trees/etc (just something small to mention in discussion).

Depth distributions

Does eye size correlate with daytime or nighttime depth distributions? Are these relationships affected by vertical migration?

Prediction: Eye size increases with daytime depth, as these animals may require the greatest sensitivity to light to cope with the exponential attenuation of light with depth.

Note: Did we have predictions about how vertical migration (which we are factoring in as the interaction effect of day and night depths) might affect eye size? I might predict that stronger migrators might have larger eyes, as they may be able to use ambient light for vision more of the time (things that stay deep at night would not have ambient light available at night). Or they might be using their eyes to follow an isolume. Just a couple ideas!

Exploratory plots (depth)

Before incorporating phylogeny, we can inspect the distribution of our data with scatter plots for absolute and relative eye size and body size vs. daytime and nighttime depth. In the below interactive plots, daytime depths are orange and nighttime depths are blue. Maximum reported depths are in black squares. Depth is not logged.

Absolute eye size with depth.

# Absolute eye size with depths ------
plot_abs_depths <- ggplot(species.phy, aes(x = eye_av, y = Day_Avg, text = genus_species)) + 
  geom_point(size = 2, alpha = 0.8, col = "orange") + 
  scale_y_reverse(name = "Depth (m)") +
  xlab("Eye diameter (mm)") + 
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  geom_point(aes(x = eye_av, y = Night_Avg), size = 2, alpha = 0.8, col = "royalblue2") +
  geom_point(aes(x = eye_av, y = Max_reported), size = 1.5, alpha = 0.8, pch = 15, col = "black")

ggplotly(plot_abs_depths)

Body size with depth.

# Body size with depth ------
plot_length_depths <- ggplot(species.phy, aes(x = length_av, y = Day_Avg, text = genus_species)) + 
  geom_point(size = 2, alpha = 0.8, col = "orange") + 
  scale_y_reverse(name = "Depth (m)") + 
  xlab("Body length (mm)") + 
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  geom_point(aes(x = length_av, y = Night_Avg), size = 2, alpha = 0.8, col = "royalblue2")+
  geom_point(aes(x = length_av, y = Max_reported), size = 1.5, alpha = 0.8, pch = 15, col = "black")

ggplotly(plot_length_depths)

Relative eye size (eye diameter/body length) with depth.

# Eye size/body size (relative eye size) with max depth -----
plot_rel_depth <- ggplot(species.phy, aes(x = eye_av/length_av, y = Day_Avg, text = genus_species)) + 
  geom_point(size = 2, alpha = 0.8, col = "orange") + 
  scale_y_reverse(name = "Depth (m)") + 
  xlab("Eye/body ratio") + 
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  geom_point(aes(x = eye_av/length_av, y = Night_Avg), size = 2, alpha = 0.8, col = "royalblue2") +
  geom_point(aes(x = eye_av/length_av, y = Max_reported), size = 1.5, alpha = 0.8, pch = 15, col = "black")

ggplotly(plot_rel_depth)

Is daytime depth correlated with maximum depth? (Presumably these are two measures of a similar thing)

#is daytime depth correlated with maximum depth?
plot_daymax <- ggplot(species.phy, aes(x = Day_Avg, y = Max_reported, text = genus_species)) + 
  geom_point(size = 2, alpha = 0.8, col = "black") + 
  scale_y_reverse(name ="Max reported depth (m)") + 
  xlab("Daytime mean depth (m)") + 
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

ggplotly(plot_daymax)

Visual inspection of the data shows that absolute eye size, body size, and relative eye size maybe could seem to increase with daytime depth, but if so it’s only driven by a couple of data points. Nighttime depth doesn’t vary much across species, but to me it looks like the species with larger body sizes might also have the larger differences between nighttime and daytime depths. Maximum reported depth shows a similar pattern to daytime depth. We next test this in an evolutionary framework.

Mean depths

Absolute eye size vs. mean day and night depths

Here, we run models of eye size vs. day x night mean depths, with lambda bounded at 0 and then again at 1.

Here is the model fit with lambda fixed at 0 (no phylogenetic signal in model residuals).

#fit model for eye diameter ~ body length (lambda = 0)
pgls_depth0 <- pgls(eye_av ~ Day_Avg * Night_Avg,
               data = species.comp, 
               lambda = 0.00001, #set lambda to 0
               bounds=list(lambda=c(0.00001,0.00001)),
               param.CI = 0.95)

#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_depth0)

par(mfrow = c(1, 1))

#main effects
anova(pgls_depth0)
## Analysis of Variance Table
## Sequential SS for pgls: lambda = 0.00, delta = 1.00, kappa = 1.00
## 
## Response: eye_av
##                   Df  Sum Sq Mean Sq F value   Pr(>F)   
## Day_Avg            1 0.60507 0.60507 10.6842 0.006720 **
## Night_Avg          1 0.00627 0.00627  0.1108 0.744988   
## Day_Avg:Night_Avg  1 0.87460 0.87460 15.4435 0.001999 **
## Residuals         12 0.67958 0.05663                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#coefficients and contrasts
summary(pgls_depth0)
## 
## Call:
## pgls(formula = eye_av ~ Day_Avg * Night_Avg, data = species.comp, 
##     lambda = 1e-05, param.CI = 0.95, bounds = list(lambda = c(1e-05, 
##         1e-05)))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.25976 -0.04300  0.05435  0.19629  0.41490 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [Fix]  : 0.000
## delta  [Fix]  : 1.000
## 
## Coefficients:
##                      Estimate  Std. Error t value Pr(>|t|)   
## (Intercept)       -3.8839e-01  3.0478e-01 -1.2744 0.226663   
## Day_Avg            1.0386e-03  2.4323e-04  4.2699 0.001088 **
## Night_Avg          6.2419e-03  1.6506e-03  3.7815 0.002617 **
## Day_Avg:Night_Avg -3.6627e-06  9.3203e-07 -3.9298 0.001999 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.238 on 12 degrees of freedom
## Multiple R-squared: 0.6862,  Adjusted R-squared: 0.6077 
## F-statistic: 8.746 on 3 and 12 DF,  p-value: 0.002394

There is a significant effect of daytime depth and a significant interaction effect of day X night depth on absolute eye size when lambda = 0.

Here is the model fit with lambda fixed at 1 (highest phylogenetic signal in model residuals).

#fit model for eye diameter ~ body length (lambda = 1)
pgls_depth1 <- pgls(eye_av ~ Day_Avg * Night_Avg,
               data = species.comp, 
               lambda = 1, #set lambda to 1
               bounds=list(lambda=c(1,1)),
               param.CI = 0.95)

#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_depth1)

par(mfrow = c(1, 1))

#main effects
anova(pgls_depth1)
## Analysis of Variance Table
## Sequential SS for pgls: lambda = 1.00, delta = 1.00, kappa = 1.00
## 
## Response: eye_av
##                   Df  Sum Sq Mean Sq F value   Pr(>F)   
## Day_Avg            1 0.05746 0.05746  1.0034 0.336259   
## Night_Avg          1 0.00178 0.00178  0.0312 0.862811   
## Day_Avg:Night_Avg  1 0.97260 0.97260 16.9856 0.001417 **
## Residuals         12 0.68712 0.05726                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#coefficients and contrasts
summary(pgls_depth1)
## 
## Call:
## pgls(formula = eye_av ~ Day_Avg * Night_Avg, data = species.comp, 
##     lambda = 1, param.CI = 0.95, bounds = list(lambda = c(1, 
##         1)))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.47134 -0.08056  0.01827  0.11141  0.30318 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [Fix]  : 1.000
## delta  [Fix]  : 1.000
## 
## Coefficients:
##                      Estimate  Std. Error t value Pr(>|t|)   
## (Intercept)       -1.0919e-01  2.7596e-01 -0.3957 0.699282   
## Day_Avg            6.6304e-04  2.0266e-04  3.2717 0.006682 **
## Night_Avg          5.8104e-03  1.4414e-03  4.0310 0.001666 **
## Day_Avg:Night_Avg -3.3074e-06  8.0249e-07 -4.1214 0.001417 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2393 on 12 degrees of freedom
## Multiple R-squared: 0.6003,  Adjusted R-squared: 0.5003 
## F-statistic: 6.007 on 3 and 12 DF,  p-value: 0.009693

When lambda = 1, only the interaction between day and night depths have a significant effect on absolute eye size.

Relative eye size vs. day and night depths.

Below, I run a PGLS model of log10(eye diameter) vs. log10(body length) with covariates of daytime depth, nighttime depth, and the interaction between daytime and nighttime depths.

#run PGLS model using the Maximum Liklihood estimate of lambda
pgls_depth <- pgls(log10(eye_av) ~ log10(length_av) + Day_Avg * Night_Avg, 
               data = species.comp, 
               lambda = "ML", bounds = list(lambda = c(1e-05, 1)),
               #uses Maximum Liklihood estimate of lambda
               param.CI = 0.95)

#evaluate model assumptions
par(mfrow = c(2,2)) #makes your plot window into 2x2 panels
plot(pgls_depth) #plot the linear model

par(mfrow = c(1,1)) #set plotting window back to one plot

The model checks here look a bit weird (residuals look bimodal; assumption is normality).

Next we can take a look at the overall significance of the main effects/interactions and sequential sums of squares. Here we can see that body length has a significant effect on eye size when depth is corrected for (expected, and we already knew that), but daytime depth, nighttime depth, and the interaction between daytime and nighttime depth do not have significant effects on relative eye size.

#anova
anova(pgls_depth)
## Analysis of Variance Table
## Sequential SS for pgls: lambda = 1.00, delta = 1.00, kappa = 1.00
## 
## Response: log10(eye_av)
##                   Df   Sum Sq  Mean Sq F value    Pr(>F)    
## log10(length_av)   1 0.253779 0.253779 74.6881 3.112e-06 ***
## Day_Avg            1 0.004141 0.004141  1.2187    0.2932    
## Night_Avg          1 0.006218 0.006218  1.8299    0.2033    
## Day_Avg:Night_Avg  1 0.000256 0.000256  0.0755    0.7886    
## Residuals         11 0.037376 0.003398                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Finally, we look at the PGLS coefficients as well significance with respect to contrasts.

#summary
summary(pgls_depth)
## 
## Call:
## pgls(formula = log10(eye_av) ~ log10(length_av) + Day_Avg * Night_Avg, 
##     data = species.comp, lambda = "ML", param.CI = 0.95, bounds = list(lambda = c(1e-05, 
##         1)))
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.075678 -0.052157  0.001252  0.039596  0.061882 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 1.000
##    lower bound : 0.000, p = 0.19099
##    upper bound : 1.000, p = 1    
##    95.0% CI   : (NA, NA)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##                      Estimate  Std. Error t value  Pr(>|t|)    
## (Intercept)       -1.3424e+00  1.7428e-01 -7.7025 9.348e-06 ***
## log10(length_av)   8.6045e-01  1.5864e-01  5.4240 0.0002089 ***
## Day_Avg            5.1702e-05  6.4331e-05  0.8037 0.4386036    
## Night_Avg         -2.4008e-04  5.7990e-04 -0.4140 0.6868326    
## Day_Avg:Night_Avg  8.7933e-08  3.2009e-07  0.2747 0.7886296    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.05829 on 11 degrees of freedom
## Multiple R-squared: 0.8761,  Adjusted R-squared: 0.8311 
## F-statistic: 19.45 on 4 and 11 DF,  p-value: 5.969e-05

The output also gives a maximum-likelihood estimation of lambda = 1.00, indicating that there is high phylogenetic signal in the residuals of the model (and that caper used a lambda of 1 to set the variance-covariance matrix in the regression). However, small sample sizes have no power to estimate the true value of lambda, so this is not very meaningful in this dataset. We should decide on a strategy for analysis (as mentioned above, may want to run all models again with lambda fixed at 0 as a comparison for the supp info so we can discuss how that affects our results).

#Likelihood plot for Pagel's lambda 
plot(pgls.profile(pgls_depth, "lambda"), 
     main = "Pagel's lambda", 
     sub = "Solid red line indicates estimate for lambda and broken red lines indcaite the 95% confidence interval")

Bounded extremes of mean depth models

With sample sizes too small to generate a robust estimate of lambda, we here again run models with the lowest and possible bounds of lambda (0 and 1), so that we can see what the range of possible parameter estimates are depending on the phylogenetic signal of the model residuals. Below, we fit these two models so that we can compare them.

Here is the model fit with lambda fixed at 0 (no phylogenetic signal in model residuals).

#fit model for eye diameter ~ body length (lambda = 0)
pgls_depth0 <- pgls(log10(eye_av) ~ log10(length_av) + Day_Avg * Night_Avg,
               data = species.comp, 
               lambda = 0.00001, #set lambda to 0
               bounds=list(lambda=c(0.00001,0.00001)),
               param.CI = 0.95)

#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_depth0)

par(mfrow = c(1, 1))

#main effects
anova(pgls_depth0)
## Analysis of Variance Table
## Sequential SS for pgls: lambda = 0.00, delta = 1.00, kappa = 1.00
## 
## Response: log10(eye_av)
##                   Df  Sum Sq Mean Sq  F value    Pr(>F)    
## log10(length_av)   1 0.35757 0.35757 113.8476 3.853e-07 ***
## Day_Avg            1 0.00801 0.00801   2.5510    0.1385    
## Night_Avg          1 0.00581 0.00581   1.8509    0.2009    
## Day_Avg:Night_Avg  1 0.00009 0.00009   0.0272    0.8719    
## Residuals         11 0.03455 0.00314                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#coefficients and contrasts
summary(pgls_depth0)
## 
## Call:
## pgls(formula = log10(eye_av) ~ log10(length_av) + Day_Avg * Night_Avg, 
##     data = species.comp, lambda = 1e-05, param.CI = 0.95, bounds = list(lambda = c(1e-05, 
##         1e-05)))
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.086086 -0.012489  0.003432  0.033701  0.096951 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [Fix]  : 0.000
## delta  [Fix]  : 1.000
## 
## Coefficients:
##                      Estimate  Std. Error t value  Pr(>|t|)    
## (Intercept)       -1.3733e+00  1.5202e-01 -9.0334 2.022e-06 ***
## log10(length_av)   8.4631e-01  1.4585e-01  5.8028 0.0001188 ***
## Day_Avg            1.1361e-04  8.0869e-05  1.4049 0.1876539    
## Night_Avg         -1.3455e-05  5.7325e-04 -0.0235 0.9816948    
## Day_Avg:Night_Avg -5.3441e-08  3.2387e-07 -0.1650 0.8719291    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.05604 on 11 degrees of freedom
## Multiple R-squared: 0.9149,  Adjusted R-squared: 0.884 
## F-statistic: 29.57 on 4 and 11 DF,  p-value: 7.848e-06

Here is the model fit with lambda fixed at 1 (highest phylogenetic signal in model residuals). *Note that when we used a ML estimate of lambda, this was the value used, though the lambda estimate was not significant.

#fit model for eye diameter ~ body length (lambda = 1)
pgls_depth1 <- pgls(log10(eye_av) ~ log10(length_av) + Day_Avg * Night_Avg,
               data = species.comp, 
               lambda = 1, #set lambda to 1
               bounds=list(lambda=c(1,1)),
               param.CI = 0.95)

#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_depth1)

par(mfrow = c(1, 1))

#main effects
anova(pgls_depth1)
## Analysis of Variance Table
## Sequential SS for pgls: lambda = 1.00, delta = 1.00, kappa = 1.00
## 
## Response: log10(eye_av)
##                   Df   Sum Sq  Mean Sq F value    Pr(>F)    
## log10(length_av)   1 0.253779 0.253779 74.6881 3.112e-06 ***
## Day_Avg            1 0.004141 0.004141  1.2187    0.2932    
## Night_Avg          1 0.006218 0.006218  1.8299    0.2033    
## Day_Avg:Night_Avg  1 0.000256 0.000256  0.0755    0.7886    
## Residuals         11 0.037376 0.003398                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#parameter estimates
summary(pgls_depth1)
## 
## Call:
## pgls(formula = log10(eye_av) ~ log10(length_av) + Day_Avg * Night_Avg, 
##     data = species.comp, lambda = 1, param.CI = 0.95, bounds = list(lambda = c(1, 
##         1)))
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.075678 -0.052157  0.001252  0.039596  0.061882 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [Fix]  : 1.000
## delta  [Fix]  : 1.000
## 
## Coefficients:
##                      Estimate  Std. Error t value  Pr(>|t|)    
## (Intercept)       -1.3424e+00  1.7428e-01 -7.7025 9.348e-06 ***
## log10(length_av)   8.6045e-01  1.5864e-01  5.4240 0.0002089 ***
## Day_Avg            5.1702e-05  6.4331e-05  0.8037 0.4386036    
## Night_Avg         -2.4008e-04  5.7990e-04 -0.4140 0.6868326    
## Day_Avg:Night_Avg  8.7933e-08  3.2009e-07  0.2747 0.7886296    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.05829 on 11 degrees of freedom
## Multiple R-squared: 0.8761,  Adjusted R-squared: 0.8311 
## F-statistic: 19.45 on 4 and 11 DF,  p-value: 5.969e-05

In both models bounding the possible true value of lambda, daytime depth, nighttime depth, and the interaction between day and night depths do not have significant effects on relative eye size in these species. This means that our finding is robust, regardless of our challenge estimating phylogenetic signal in our model residuals.

Maximum depth

While we found no association between eye size and mean daytime or nighttime depths, we can also test eye size vs. the maximum reported depth for each species.

Absolute eye size vs. max depth

Here is a model of eye size ~ maximum reported depth fit with lambda fixed at 0 (no phylogenetic signal in model residuals).

#fit model for eye diameter ~ body length (lambda = 0)
pgls_maxdepth0 <- pgls(eye_av ~ Max_reported,
               data = species.comp, 
               lambda = 0.00001, #set lambda to 0
               bounds=list(lambda=c(0.00001,0.00001)),
               param.CI = 0.95)

#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_maxdepth0)

par(mfrow = c(1, 1))

#main effects
anova(pgls_maxdepth0)
## Analysis of Variance Table
## Sequential SS for pgls: lambda = 0.00, delta = 1.00, kappa = 1.00
## 
## Response: eye_av
##              Df Sum Sq Mean Sq F value   Pr(>F)   
## Max_reported  1 1.0314 1.03141  12.732 0.003087 **
## Residuals    14 1.1341 0.08101                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#coefficients and contrasts
summary(pgls_maxdepth0)
## 
## Call:
## pgls(formula = eye_av ~ Max_reported, data = species.comp, lambda = 1e-05, 
##     param.CI = 0.95, bounds = list(lambda = c(1e-05, 1e-05)))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.35034 -0.20774 -0.03950  0.05696  0.58971 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [Fix]  : 0.000
## delta  [Fix]  : 1.000
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  0.26375005 0.23390646  1.1276 0.278449   
## Max_reported 0.00049307 0.00013818  3.5682 0.003087 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2846 on 14 degrees of freedom
## Multiple R-squared: 0.4763,  Adjusted R-squared: 0.4389 
## F-statistic: 12.73 on 1 and 14 DF,  p-value: 0.003087

Maximum reported depth has a significant effect on absolute eye size when lambda = 0.

Here is the model fit with lambda fixed at 1 (highest phylogenetic signal in model residuals).

#fit model for eye diameter ~ body length (lambda = 1)
pgls_maxdepth1 <- pgls(eye_av ~ Max_reported,
               data = species.comp, 
               lambda = 1, #set lambda to 1
               bounds=list(lambda=c(1,1)),
               param.CI = 0.95)

#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_maxdepth1)

par(mfrow = c(1, 1))

#main effects
anova(pgls_maxdepth1)
## Analysis of Variance Table
## Sequential SS for pgls: lambda = 1.00, delta = 1.00, kappa = 1.00
## 
## Response: eye_av
##              Df  Sum Sq Mean Sq F value  Pr(>F)  
## Max_reported  1 0.45511 0.45511  5.0414 0.04142 *
## Residuals    14 1.26385 0.09028                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#parameter estimates
summary(pgls_maxdepth1)
## 
## Call:
## pgls(formula = eye_av ~ Max_reported, data = species.comp, lambda = 1, 
##     param.CI = 0.95, bounds = list(lambda = c(1, 1)))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.66883 -0.28878 -0.07065  0.02554  0.37075 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [Fix]  : 1.000
## delta  [Fix]  : 1.000
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  0.42595485 0.27755652  1.5347  0.14715  
## Max_reported 0.00036789 0.00016385  2.2453  0.04142 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3005 on 14 degrees of freedom
## Multiple R-squared: 0.2648,  Adjusted R-squared: 0.2122 
## F-statistic: 5.041 on 1 and 14 DF,  p-value: 0.04142

Maximum reported depth has a significant effect on absolute eye size when lambda = 1.

In both models bounding the possible true value of lambda, maximum depth has a significant effect on absolute eye size.

Relative eye size vs. max depth

Below, I run a PGLS model of log10(eye diameter) vs. log10(body length) with covariate of maximum depth.

#run PGLS model using the Maximum Liklihood estimate of lambda
pgls_maxdepth <- pgls(log10(eye_av) ~ log10(length_av)  + Max_reported, 
               data = species.comp, 
               lambda = "ML", bounds = list(lambda = c(1e-05, 1)),
               #uses Maximum Liklihood estimate of lambda
               param.CI = 0.95)

#evaluate model assumptions
par(mfrow = c(2,2)) #makes your plot window into 2x2 panels
plot(pgls_maxdepth) #plot the linear model

par(mfrow = c(1,1)) #set plotting window back to one plot

The model checks here look ok.

Next we can take a look at the overall significance of the main effects/interactions and sequential sums of squares. Here we can see that body length has a significant effect on eye size when maximum reported depth is corrected for (expected, and we already knew that), but that maximum depth does not have a significant effect on eye size.

#anova
anova(pgls_maxdepth)
## Analysis of Variance Table
## Sequential SS for pgls: lambda = 0.89, delta = 1.00, kappa = 1.00
## 
## Response: log10(eye_av)
##                  Df   Sum Sq  Mean Sq F value    Pr(>F)    
## log10(length_av)  1 0.257318 0.257318  82.669 5.364e-07 ***
## Max_reported      1 0.006247 0.006247   2.007    0.1801    
## Residuals        13 0.040464 0.003113                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Finally, we look at the PGLS coefficients as well significance with respect to contrasts.

#summary
summary(pgls_maxdepth)
## 
## Call:
## pgls(formula = log10(eye_av) ~ log10(length_av) + Max_reported, 
##     data = species.comp, lambda = "ML", param.CI = 0.95, bounds = list(lambda = c(1e-05, 
##         1)))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.05638 -0.02305  0.01184  0.06568  0.08418 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 0.886
##    lower bound : 0.000, p = 0.28682
##    upper bound : 1.000, p = 0.80796
##    95.0% CI   : (NA, NA)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##                     Estimate  Std. Error t value  Pr(>|t|)    
## (Intercept)      -1.2799e+00  1.4323e-01 -8.9360 6.528e-07 ***
## log10(length_av)  7.6957e-01  1.0295e-01  7.4751 4.663e-06 ***
## Max_reported      4.8325e-05  3.4111e-05  1.4167    0.1801    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.05579 on 13 degrees of freedom
## Multiple R-squared: 0.8669,  Adjusted R-squared: 0.8464 
## F-statistic: 42.34 on 2 and 13 DF,  p-value: 2.028e-06

The output also gives a maximum-likelihood estimation of lambda = 0.89, indicating that there is high phylogenetic signal in the residuals of the model (and that caper used a lambda of 0.89 to set the variance-covariance matrix in the regression). However, small sample sizes have no power to estimate the true value of lambda, so this is not very meaningful in this dataset. We additionally run models with lambda fixed at 0 and at 1 as comparisons for the supp info so we can discuss how that affects our results.

#Likelihood plot for Pagel's lambda 
plot(pgls.profile(pgls_maxdepth, "lambda"), 
     main = "Pagel's lambda", 
     sub = "Solid red line indicates estimate for lambda and broken red lines indcaite the 95% confidence interval")

Bounded extremes of maximum reported depth models

With sample sizes too small to generate a robust estimate of lambda, we here run models with the lowest and possible bounds of lambda (0 and 1), so that we can see what the range of possible results are depending on the phylogenetic signal of the model residuals. Below, we fit these two models so that we can compare them.

Here is the model fit with lambda fixed at 0 (no phylogenetic signal in model residuals).

#fit model for eye diameter ~ body length (lambda = 0)
pgls_maxdepth0 <- pgls(log10(eye_av) ~ log10(length_av) + Max_reported,
               data = species.comp, 
               lambda = 0.00001, #set lambda to 0
               bounds=list(lambda=c(0.00001,0.00001)),
               param.CI = 0.95)

#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_maxdepth0)

par(mfrow = c(1, 1))

#main effects
anova(pgls_maxdepth0)
## Analysis of Variance Table
## Sequential SS for pgls: lambda = 0.00, delta = 1.00, kappa = 1.00
## 
## Response: log10(eye_av)
##                  Df  Sum Sq Mean Sq  F value    Pr(>F)    
## log10(length_av)  1 0.35757 0.35757 123.7525 5.157e-08 ***
## Max_reported      1 0.01090 0.01090   3.7715   0.07412 .  
## Residuals        13 0.03756 0.00289                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#parameter estimates
summary(pgls_maxdepth0)
## 
## Call:
## pgls(formula = log10(eye_av) ~ log10(length_av) + Max_reported, 
##     data = species.comp, lambda = 1e-05, param.CI = 0.95, bounds = list(lambda = c(1e-05, 
##         1e-05)))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.09807 -0.03771  0.01952  0.04081  0.07693 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [Fix]  : 0.000
## delta  [Fix]  : 1.000
## 
## Coefficients:
##                     Estimate  Std. Error t value  Pr(>|t|)    
## (Intercept)      -1.3396e+00  1.3743e-01 -9.7469 2.421e-07 ***
## log10(length_av)  8.0134e-01  1.0607e-01  7.5547 4.162e-06 ***
## Max_reported      6.4411e-05  3.3167e-05  1.9420   0.07412 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.05375 on 13 degrees of freedom
## Multiple R-squared: 0.9075,  Adjusted R-squared: 0.8933 
## F-statistic: 63.76 on 2 and 13 DF,  p-value: 1.907e-07

When lambda = 0, maximum depth has no effect on relative eye size.

Here is the model fit with lambda fixed at 1 (highest phylogenetic signal in model residuals).

#fit model for eye diameter ~ body length (lambda = 1)
pgls_maxdepth1 <- pgls(log10(eye_av) ~ log10(length_av) + Max_reported,
               data = species.comp, 
               lambda = 1, #set lambda to 1
               bounds=list(lambda=c(1,1)),
               param.CI = 0.95)

#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_maxdepth1)

par(mfrow = c(1, 1))

#main effects
anova(pgls_maxdepth1)
## Analysis of Variance Table
## Sequential SS for pgls: lambda = 1.00, delta = 1.00, kappa = 1.00
## 
## Response: log10(eye_av)
##                  Df   Sum Sq  Mean Sq F value    Pr(>F)    
## log10(length_av)  1 0.253779 0.253779 78.0290 7.436e-07 ***
## Max_reported      1 0.005711 0.005711  1.7559     0.208    
## Residuals        13 0.042281 0.003252                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#parameter estimates
summary(pgls_maxdepth1)
## 
## Call:
## pgls(formula = log10(eye_av) ~ log10(length_av) + Max_reported, 
##     data = species.comp, lambda = 1, param.CI = 0.95, bounds = list(lambda = c(1, 
##         1)))
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.072284 -0.035744  0.008099  0.045397  0.078221 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [Fix]  : 1.000
## delta  [Fix]  : 1.000
## 
## Coefficients:
##                     Estimate  Std. Error t value  Pr(>|t|)    
## (Intercept)      -1.2672e+00  1.4497e-01 -8.7411 8.369e-07 ***
## log10(length_av)  7.6322e-01  1.0345e-01  7.3773 5.367e-06 ***
## Max_reported      4.5769e-05  3.4539e-05  1.3251     0.208    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.05703 on 13 degrees of freedom
## Multiple R-squared: 0.8599,  Adjusted R-squared: 0.8383 
## F-statistic: 39.89 on 2 and 13 DF,  p-value: 2.832e-06

When lambda = 1, maximum depth has no effect on relative eye size.

In both models bounding the possible true value of lambda, maximum depth does not have significant effects on relative eye size in these species.

Depth range traversed in DVM (daytime mean - nighttime mean)

While we found no association between eye size and mean daytime or nighttime depths or the interaction between them, I also wanted to check whether there was an association between eye size and the distance traversed during DVM.

Absolute eye size vs. depths traversed

Here is the model of eye size vs. depth range fit with lambda fixed at 0 (no phylogenetic signal in model residuals).

#fit model for eye diameter ~ body length (lambda = 0)
pgls_depthrange0 <- pgls(eye_av ~ Range,
               data = species.comp, 
               lambda = 0.00001, #set lambda to 0
               bounds=list(lambda=c(0.00001,0.00001)),
               param.CI = 0.95)

#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_depthrange0)

par(mfrow = c(1, 1))

#main effects
anova(pgls_depthrange0)
## Analysis of Variance Table
## Sequential SS for pgls: lambda = 0.00, delta = 1.00, kappa = 1.00
## 
## Response: eye_av
##           Df Sum Sq Mean Sq F value  Pr(>F)  
## Range      1 0.4582 0.45820  3.7572 0.07301 .
## Residuals 14 1.7073 0.12195                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#parameter estimates
summary(pgls_depthrange0)
## 
## Call:
## pgls(formula = eye_av ~ Range, data = species.comp, lambda = 1e-05, 
##     param.CI = 0.95, bounds = list(lambda = c(1e-05, 1e-05)))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.62839 -0.26817  0.04641  0.22210  0.56496 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [Fix]  : 0.000
## delta  [Fix]  : 1.000
## 
## Coefficients:
##               Estimate Std. Error t value  Pr(>|t|)    
## (Intercept) 0.78225824 0.16726887  4.6767 0.0003566 ***
## Range       0.00055569 0.00028668  1.9384 0.0730117 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3492 on 14 degrees of freedom
## Multiple R-squared: 0.2116,  Adjusted R-squared: 0.1553 
## F-statistic: 3.757 on 1 and 14 DF,  p-value: 0.07301

When lambda = 0, depth range has no significant effect on absolute eye size.

Here is the model fit with lambda fixed at 1 (highest phylogenetic signal in model residuals).

#fit model for eye diameter ~ body length (lambda = 1)
pgls_depthrange1 <- pgls(eye_av ~ Range,
               data = species.comp, 
               lambda = 1, #set lambda to 1
               bounds=list(lambda=c(1,1)),
               param.CI = 0.95)

#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_depthrange1)

par(mfrow = c(1, 1))

#main effects
anova(pgls_depthrange1)
## Analysis of Variance Table
## Sequential SS for pgls: lambda = 1.00, delta = 1.00, kappa = 1.00
## 
## Response: eye_av
##           Df  Sum Sq  Mean Sq F value Pr(>F)
## Range      1 0.04454 0.044539  0.3724 0.5515
## Residuals 14 1.67442 0.119602
#parameter estimates
summary(pgls_depthrange1)
## 
## Call:
## pgls(formula = eye_av ~ Range, data = species.comp, lambda = 1, 
##     param.CI = 0.95, bounds = list(lambda = c(1, 1)))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.64588 -0.44759 -0.13224  0.03535  0.27527 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [Fix]  : 1.000
## delta  [Fix]  : 1.000
## 
## Coefficients:
##               Estimate Std. Error t value  Pr(>|t|)    
## (Intercept) 0.93155235 0.16629239  5.6019 6.524e-05 ***
## Range       0.00013888 0.00022758  0.6102    0.5515    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3458 on 14 degrees of freedom
## Multiple R-squared: 0.02591, Adjusted R-squared: -0.04367 
## F-statistic: 0.3724 on 1 and 14 DF,  p-value: 0.5515

When lambda = 0, depth range has no significant effect on absolute eye size. Both models agree there is no effect of depth traversed during migration on absolute eye size.

Relative eye size vs. depth range traversed

Below, I run a PGLS model of log10(eye diameter) vs. log10(body length) with depth range traversed during DVM as a covariate.

#run PGLS model using the Maximum Liklihood estimate of lambda
pgls_depthrange <- pgls(log10(eye_av) ~ log10(length_av) + Range, 
               data = species.comp, 
               lambda = "ML", bounds = list(lambda = c(1e-05, 1)),
               #uses Maximum Liklihood estimate of lambda
               param.CI = 0.95)

#evaluate model assumptions
par(mfrow = c(2,2)) #makes your plot window into 2x2 panels
plot(pgls_depthrange) #plot the linear model

par(mfrow = c(1,1)) #set plotting window back to one plot

The model checks here look a bit wonky.

Next we can take a look at the overall significance of the main effects and sequential sums of squares. Here we can see that body length has a significant effect on eye size when depth range is corrected for (expected, and we already knew that), but that depth range does not have a significant effect on eye size.

#anova
anova(pgls_depthrange)
## Analysis of Variance Table
## Sequential SS for pgls: lambda = 1.00, delta = 1.00, kappa = 1.00
## 
## Response: log10(eye_av)
##                  Df   Sum Sq  Mean Sq F value    Pr(>F)    
## log10(length_av)  1 0.253779 0.253779 86.7481 4.078e-07 ***
## Range             1 0.009961 0.009961  3.4048    0.0879 .  
## Residuals        13 0.038031 0.002925                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Finally, we look at the PGLS coefficients as well significance with respect to contrasts.

#summary
summary(pgls_depthrange)
## 
## Call:
## pgls(formula = log10(eye_av) ~ log10(length_av) + Range, data = species.comp, 
##     lambda = "ML", param.CI = 0.95, bounds = list(lambda = c(1e-05, 
##         1)))
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.069275 -0.051016  0.008185  0.041543  0.058192 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 1.000
##    lower bound : 0.000, p = 0.22392
##    upper bound : 1.000, p = 1    
##    95.0% CI   : (NA, NA)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##                     Estimate  Std. Error t value  Pr(>|t|)    
## (Intercept)      -1.3117e+00  1.3728e-01 -9.5547 3.044e-07 ***
## log10(length_av)  8.1984e-01  8.8362e-02  9.2782 4.261e-07 ***
## Range             6.5688e-05  3.5599e-05  1.8452    0.0879 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.05409 on 13 degrees of freedom
## Multiple R-squared: 0.874,   Adjusted R-squared: 0.8546 
## F-statistic: 45.08 on 2 and 13 DF,  p-value: 1.422e-06

The output also gives a maximum-likelihood estimation of lambda = 1, indicating that there is high phylogenetic signal in the residuals of the model (and that caper used a lambda of 1 to set the variance-covariance matrix in the regression). However, small sample sizes have no power to estimate the true value of lambda, so this is not very meaningful in this dataset. We additionally run models with lambda fixed at 0 and at 1 as comparisons for the supp info so we can discuss how that affects our results.

#Likelihood plot for Pagel's lambda 
plot(pgls.profile(pgls_depthrange, "lambda"), 
     main = "Pagel's lambda", 
     sub = "Solid red line indicates estimate for lambda and broken red lines indcaite the 95% confidence interval")

Bounded extremes of depth range models

With sample sizes too small to generate a robust estimate of lambda, we here run models with the lowest and highest possible bounds of lambda (0 and 1), so that we can see what the range of possible results are depending on the phylogenetic signal of the model residuals. Below, we fit these two models so that we can compare them.

Here is the model fit with lambda fixed at 0 (no phylogenetic signal in model residuals).

#fit model for eye diameter ~ body length (lambda = 0)
pgls_depthrange0 <- pgls(log10(eye_av) ~ log10(length_av) + Range,
               data = species.comp, 
               lambda = 0.00001, #set lambda to 0
               bounds=list(lambda=c(0.00001,0.00001)),
               param.CI = 0.95)

#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_depthrange0)

par(mfrow = c(1, 1))

#main effects
anova(pgls_depthrange0)
## Analysis of Variance Table
## Sequential SS for pgls: lambda = 0.00, delta = 1.00, kappa = 1.00
## 
## Response: log10(eye_av)
##                  Df  Sum Sq Mean Sq F value    Pr(>F)    
## log10(length_av)  1 0.35757 0.35757 134.153 3.192e-08 ***
## Range             1 0.01381 0.01381   5.181    0.0404 *  
## Residuals        13 0.03465 0.00267                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#parameter estimates
summary(pgls_depthrange0)
## 
## Call:
## pgls(formula = log10(eye_av) ~ log10(length_av) + Range, data = species.comp, 
##     lambda = 1e-05, param.CI = 0.95, bounds = list(lambda = c(1e-05, 
##         1e-05)))
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.088931 -0.012322  0.002289  0.031130  0.099003 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [Fix]  : 0.000
## delta  [Fix]  : 1.000
## 
## Coefficients:
##                     Estimate  Std. Error  t value  Pr(>|t|)    
## (Intercept)      -1.3794e+00  1.2604e-01 -10.9445 6.253e-08 ***
## log10(length_av)  8.6156e-01  8.5383e-02  10.0905 1.621e-07 ***
## Range             1.0275e-04  4.5143e-05   2.2762    0.0404 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.05163 on 13 degrees of freedom
## Multiple R-squared: 0.9147,  Adjusted R-squared: 0.9015 
## F-statistic: 69.67 on 2 and 13 DF,  p-value: 1.128e-07

When lambda = 0, depth range traversed has a significant effect on relative eye size.

Here is the model fit with lambda fixed at 1 (highest phylogenetic signal in model residuals).

#fit model for eye diameter ~ body length (lambda = 1)
pgls_depthrange1 <- pgls(log10(eye_av) ~ log10(length_av) + Range,
               data = species.comp, 
               lambda = 1, #set lambda to 1
               bounds=list(lambda=c(1,1)),
               param.CI = 0.95)

#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_depthrange1)

par(mfrow = c(1, 1))

#main effects
anova(pgls_depthrange1)
## Analysis of Variance Table
## Sequential SS for pgls: lambda = 1.00, delta = 1.00, kappa = 1.00
## 
## Response: log10(eye_av)
##                  Df   Sum Sq  Mean Sq F value    Pr(>F)    
## log10(length_av)  1 0.253779 0.253779 86.7481 4.078e-07 ***
## Range             1 0.009961 0.009961  3.4048    0.0879 .  
## Residuals        13 0.038031 0.002925                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#parameter estimates
summary(pgls_depthrange1)
## 
## Call:
## pgls(formula = log10(eye_av) ~ log10(length_av) + Range, data = species.comp, 
##     lambda = 1, param.CI = 0.95, bounds = list(lambda = c(1, 
##         1)))
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.069275 -0.051016  0.008185  0.041543  0.058192 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [Fix]  : 1.000
## delta  [Fix]  : 1.000
## 
## Coefficients:
##                     Estimate  Std. Error t value  Pr(>|t|)    
## (Intercept)      -1.3117e+00  1.3728e-01 -9.5547 3.044e-07 ***
## log10(length_av)  8.1984e-01  8.8362e-02  9.2782 4.261e-07 ***
## Range             6.5688e-05  3.5599e-05  1.8452    0.0879 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.05409 on 13 degrees of freedom
## Multiple R-squared: 0.874,   Adjusted R-squared: 0.8546 
## F-statistic: 45.08 on 2 and 13 DF,  p-value: 1.422e-06

In the model with lambda set to 0, depth range does not have a significant effect on relative eye size. Thus, the answer depends on the true value of lambda (which we can’t get at with our small sample size).

Let’s just see what happens with a lambda of 0.5 for better understanding.

#fit model for eye diameter ~ body length (lambda = 1)
pgls_depthrange05 <- pgls(log10(eye_av) ~ log10(length_av) + Range,
               data = species.comp, 
               lambda = 0.5, #set lambda to 0.5
               param.CI = 0.95)

#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_depthrange1)

par(mfrow = c(1, 1))

#main effects
anova(pgls_depthrange05)
## Analysis of Variance Table
## Sequential SS for pgls: lambda = 0.50, delta = 1.00, kappa = 1.00
## 
## Response: log10(eye_av)
##                  Df   Sum Sq  Mean Sq  F value    Pr(>F)    
## log10(length_av)  1 0.280831 0.280831 103.4663 1.476e-07 ***
## Range             1 0.010073 0.010073   3.7112    0.0762 .  
## Residuals        13 0.035285 0.002714                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#parameter estimates
summary(pgls_depthrange05)
## 
## Call:
## pgls(formula = log10(eye_av) ~ log10(length_av) + Range, data = species.comp, 
##     lambda = 0.5, param.CI = 0.95)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.06357 -0.02033  0.03424  0.05163  0.08144 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [Fix]  : 0.500
## delta  [Fix]  : 1.000
## 
## Coefficients:
##                     Estimate  Std. Error  t value  Pr(>|t|)    
## (Intercept)      -1.3526e+00  1.3250e-01 -10.2085 1.415e-07 ***
## log10(length_av)  8.4502e-01  8.7069e-02   9.7052 2.543e-07 ***
## Range             8.2101e-05  4.2618e-05   1.9264    0.0762 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0521 on 13 degrees of freedom
## Multiple R-squared: 0.8918,  Adjusted R-squared: 0.8752 
## F-statistic: 53.59 on 2 and 13 DF,  p-value: 5.27e-07

At lambda = 0.5 depth range does not have a significant effect on relative eye size.

Habitat

Preliminary plots (habitat)

Before incorporating phylogeny, we can inspect the distribution of our data with boxplots for absolute and relative eye size be habitat.

# Phylogeny with eye size and investment and ecology -------

#color vector for photophores
col_eco <- c("no" = "#7FDBFF",
            "yes" = "#FF851B")

#resort trait dataset to the order of tree tip labels
species.phy <- species.phy[tree.plot$tip.label, ] 

#plot tree with eye investment bars
plotTree.wBars(tree.plot, releye, 
               scale = 1.5, 
               tip.labels = TRUE,
               col = unname(col_eco[as.factor(species.phy$Benthopelagic)]),
               offset = 1,
               ftype = "bi",
               fsize = 0.7)

#add tip labels for absolute eye size
tiplabels(cex = 1.5*aveye,
          pch = 19, #shape of labels
          col = "mediumpurple",
          offset = 0)

#add legend for photophores
legend(x = "bottomleft", legend = c("pelagic", "benthopelagic"), pch = 22, pt.cex= 2, pt.bg = col_eco, cex = 1, bty = "n", horiz = F)

#Boxplots for absolute eye size across photophores -----
plot_abs_eco <- ggplot(data = species.phy, aes(x = reorder(Benthopelagic, eye_av, FUN = mean), y = eye_av)) + 
  geom_boxplot(notch = FALSE, outlier.shape = NA, alpha = 0.6) + #controls boxes
  stat_summary(fun.y = mean, colour = "black", geom = "point", shape = 18, size = 3, show_guide = FALSE) + #controls what stats shown
 geom_jitter(shape = 19, size = 1.3, alpha = 0.5, position = position_jitter(0.15)) +
  theme(text = element_text(size=14), panel.background = element_blank(), axis.line = element_line(colour = "black")) + #controls background
  xlab("Benthopelagic") +
  ylab("Eye diameter (mm)") +
  ggtitle("Absolute eye size")

ggplotly(plot_abs_eco)
#Boxplots for relative eye investment across photophores -----
plot_rel_eco <- ggplot(data = species.phy, aes(x = reorder(Benthopelagic, pglsres, FUN = mean), y = pglsres)) + 
  geom_boxplot(notch = FALSE, outlier.shape = NA, alpha = 0.6) + #controls boxes
  stat_summary(fun.y = mean, colour = "black", geom = "point", shape = 18, size = 3, show_guide = FALSE) + #controls what stats shown
 geom_jitter(shape = 19, size = 1.3, alpha = 0.5, position = position_jitter(0.15)) +
  theme(text = element_text(size=14), panel.background = element_blank(), axis.line = element_line(colour = "black")) + #controls background
  xlab("Benthopelagic") +
  ylab("Eye ~ body residuals") +
  ggtitle("Relative eye investment")

ggplotly(plot_rel_eco) 
#Export figure------
pdf("../Figures/phylo_eco.pdf", width = 6, height = 4)

#plot tree with eye investment bars
plotTree.wBars(tree.plot, releye, 
               scale = 1.5, 
               tip.labels = TRUE,
               col = unname(col_eco[as.factor(species.phy$Benthopelagic)]),
               offset = 1,
               ftype = "bi",
               fsize = 0.7)


#add tip labels for absolute eye size
tiplabels(cex = 1.5*aveye,
          pch = 19, #shape of labels
          col = "mediumpurple",
          offset = 0)

#add legend for ecology
legend(x = "bottomleft", legend = c("pelagic", "benthopelagic"), pch = 22, pt.cex= 2, pt.bg = col_eco, cex = 1, bty = "n", horiz = F)

dev.off()

It looks like benthopelagic species have big-ish eyes, but it’s small sample size so hard to say. Distribution overlaps with pelagic species. We next test this in an evolutionary framework.

Absolute eye size vs. habitat

We test for effects of habitat on absolute eye size using a PGLS of eye size ~ habitat.

Here is a model fit with lambda fixed at 0 (no phylogenetic signal in model residuals).

#fit model for eye diameter ~ body length (lambda = 0)
pgls_eco0 <- pgls(eye_av ~ as.factor(Benthopelagic), 
               data = species.comp, 
               lambda = 0.00001, #set lambda to 0
               bounds=list(lambda=c(0.00001,0.00001)),
               param.CI = 0.95)

#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_eco0)

par(mfrow = c(1, 1))

#main effects
anova(pgls_eco0)
## Analysis of Variance Table
## Sequential SS for pgls: lambda = 0.00, delta = 1.00, kappa = 1.00
## 
## Response: eye_av
##                          Df  Sum Sq Mean Sq F value Pr(>F)
## as.factor(Benthopelagic)  1 0.11959 0.11959  0.8183  0.381
## Residuals                14 2.04594 0.14614
#parameter estimates
summary(pgls_eco0)
## 
## Call:
## pgls(formula = eye_av ~ as.factor(Benthopelagic), data = species.comp, 
##     lambda = 1e-05, param.CI = 0.95, bounds = list(lambda = c(1e-05, 
##         1e-05)))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.59327 -0.29244  0.02655  0.26145  0.58982 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [Fix]  : 0.000
## delta  [Fix]  : 1.000
## 
## Coefficients:
##                             Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)                  1.01729    0.10603  9.5947 1.556e-07 ***
## as.factor(Benthopelagic)yes  0.22150    0.24486  0.9046     0.381    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3823 on 14 degrees of freedom
## Multiple R-squared: 0.05522, Adjusted R-squared: -0.01226 
## F-statistic: 0.8183 on 1 and 14 DF,  p-value: 0.381

When lambda = 0, habitat has no effect on absolute eye size.

Here is the model fit with lambda fixed at 1 (highest phylogenetic signal in model residuals).

#fit model for eye diameter ~ body length (lambda = 1)
pgls_eco1 <- pgls(eye_av ~ as.factor(Benthopelagic), 
               data = species.comp, 
               lambda = 1, #set lambda to 1
               bounds=list(lambda=c(1,1)),
               param.CI = 0.95)

#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_eco1)

par(mfrow = c(1, 1))

#main effects
anova(pgls_eco1)
## Analysis of Variance Table
## Sequential SS for pgls: lambda = 1.00, delta = 1.00, kappa = 1.00
## 
## Response: eye_av
##                          Df  Sum Sq  Mean Sq F value Pr(>F)
## as.factor(Benthopelagic)  1 0.00128 0.001281  0.0104 0.9201
## Residuals                14 1.71768 0.122692
#parameter estimates
summary(pgls_eco1)
## 
## Call:
## pgls(formula = eye_av ~ as.factor(Benthopelagic), data = species.comp, 
##     lambda = 1, param.CI = 0.95, bounds = list(lambda = c(1, 
##         1)))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.71189 -0.44229 -0.13957  0.03058  0.32764 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [Fix]  : 1.000
## delta  [Fix]  : 1.000
## 
## Coefficients:
##                             Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)                 0.988406   0.138957  7.1130 5.229e-06 ***
## as.factor(Benthopelagic)yes 0.021808   0.213436  0.1022    0.9201    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3503 on 14 degrees of freedom
## Multiple R-squared: 0.0007451,   Adjusted R-squared: -0.07063 
## F-statistic: 0.01044 on 1 and 14 DF,  p-value: 0.9201

Habitat (benthopelagic/pelagic) does not have a significant effect on absolute eye size when lambda is set to 0 or 1. Thus, our finding is robust whatever the true value of lambda is. However, with n = 3 for benthopelagic species, we should note that our power here is very low.

Relative eye size vs. habitat

Below, I run a PGLS model of log10(eye diameter) vs. log10(body length) with habitat (benthopelagic vs. pelagic) as a covariate.

#run PGLS modelusing the Maximum Liklihood estimate of lambda
pgls_eco <- pgls(log10(eye_av) ~ log10(length_av) + Benthopelagic, 
               data = species.comp, 
               lambda = "ML",
               #uses Maximum Liklihood estimate of lambda
               param.CI = 0.95)

#evaluate model assumptions
par(mfrow = c(2,2)) #makes your plot window into 2x2 panels
plot(pgls_eco) #plot the linear model

par(mfrow = c(1,1)) #set plotting window back to one plot

Next we can take a look at the overall significance of the main, 2-way, 3-way, etc. interactions and sequential sums of squares. Here we can see that body length has a significant effect on eye size (expected, we already knew that), but ecology (benthopelagic or not) has no effect on eye size (p = 0.57).

#anova
anova(pgls_eco)
## Analysis of Variance Table
## Sequential SS for pgls: lambda = 1.00, delta = 1.00, kappa = 1.00
## 
## Response: log10(eye_av)
##                  Df   Sum Sq  Mean Sq F value    Pr(>F)    
## log10(length_av)  1 0.253779 0.253779 70.4597 1.316e-06 ***
## Benthopelagic     1 0.001169 0.001169  0.3245    0.5786    
## Residuals        13 0.046823 0.003602                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Finally, we look at the PGLS coefficients as well as significance with respect to contrasts.

#summary
summary(pgls_eco)
## 
## Call:
## pgls(formula = log10(eye_av) ~ log10(length_av) + Benthopelagic, 
##     data = species.comp, lambda = "ML", param.CI = 0.95)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.092420 -0.049981  0.009349  0.043498  0.080494 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 1.000
##    lower bound : 0.000, p = 0.11794
##    upper bound : 1.000, p = 1    
##    95.0% CI   : (NA, NA)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##                   Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)      -1.287737   0.151661 -8.4909 1.158e-06 ***
## log10(length_av)  0.820905   0.098088  8.3691 1.360e-06 ***
## Benthopelagicyes  0.020845   0.036592  0.5697    0.5786    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.06001 on 13 degrees of freedom
## Multiple R-squared: 0.8448,  Adjusted R-squared: 0.821 
## F-statistic: 35.39 on 2 and 13 DF,  p-value: 5.496e-06

The output also includes the maximum-likelihood estimation of lambda, which represents the phylogenetic signal for PGLS model residuals. Here, the estimate for lambda is 1.00, which means that there is the highest phylogenetic signal. However, as mentioned above, when you have <20 species the power to estimate lambda is not there, so below we will run models with the min and max possible values of lambda to see how that affects our results.

#Likelihood plot for Pagel's lambda 
plot(pgls.profile(pgls_eco, "lambda"), 
     main = "Pagel's lambda", 
     sub = "Solid red line indicates estimate for lambda and broken red lines indcaite the 95% confidence interval")

Bounded extremes of habitat models

With sample sizes too small to generate a robust estimate of lambda, one option is to run models with the lowest and possible bounds of lambda (0 and 1), so that we can see what the range of possible parameter estimates are depending on the phylogenetic signal of the model residuals. Below, we fit these two models so that we can compare them.

Here is the model fit with lambda fixed at 0 (no phylogenetic signal in model residuals).

#fit model for eye diameter ~ body length (lambda = 0)
pgls_eco0 <- pgls(log10(eye_av) ~ log10(length_av) + Benthopelagic, 
               data = species.comp, 
               lambda = 0.00001, #set lambda to 0
               bounds=list(lambda=c(0.00001,0.00001)),
               param.CI = 0.95)

#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_eco0)

par(mfrow = c(1, 1))

#main effects
anova(pgls_eco0)
## Analysis of Variance Table
## Sequential SS for pgls: lambda = 0.00, delta = 1.00, kappa = 1.00
## 
## Response: log10(eye_av)
##                  Df  Sum Sq Mean Sq  F value    Pr(>F)    
## log10(length_av)  1 0.35757 0.35757 102.5837 1.551e-07 ***
## Benthopelagic     1 0.00315 0.00315   0.9026    0.3594    
## Residuals        13 0.04531 0.00349                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#parameter estimates
summary(pgls_eco0)
## 
## Call:
## pgls(formula = log10(eye_av) ~ log10(length_av) + Benthopelagic, 
##     data = species.comp, lambda = 1e-05, param.CI = 0.95, bounds = list(lambda = c(1e-05, 
##         1e-05)))
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.102084 -0.026881  0.005142  0.048346  0.095391 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [Fix]  : 0.000
## delta  [Fix]  : 1.000
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -1.409336   0.143626 -9.8125 2.24e-07 ***
## log10(length_av)  0.909764   0.093762  9.7029 2.55e-07 ***
## Benthopelagicyes  0.036746   0.038678  0.9501   0.3594    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.05904 on 13 degrees of freedom
## Multiple R-squared: 0.8884,  Adjusted R-squared: 0.8712 
## F-statistic: 51.74 on 2 and 13 DF,  p-value: 6.454e-07

With lambda = 0, habitat has no effect on relative eye size.

Here is the model fit with lambda fixed at 1 (highest phylogenetic signal in model residuals). *Note that when we used a ML estimate of lambda, this was the value used, though the lambda estimate was not significant.

#fit model for eye diameter ~ body length (lambda = 1)
pgls_eco1 <- pgls(log10(eye_av) ~ log10(length_av) + Benthopelagic, 
               data = species.comp, 
               lambda = 1, #set lambda to 1
               bounds=list(lambda=c(1,1)),
               param.CI = 0.95)

#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_eco1)

par(mfrow = c(1, 1))

#main effects
anova(pgls_eco1)
## Analysis of Variance Table
## Sequential SS for pgls: lambda = 1.00, delta = 1.00, kappa = 1.00
## 
## Response: log10(eye_av)
##                  Df   Sum Sq  Mean Sq F value    Pr(>F)    
## log10(length_av)  1 0.253779 0.253779 70.4597 1.316e-06 ***
## Benthopelagic     1 0.001169 0.001169  0.3245    0.5786    
## Residuals        13 0.046823 0.003602                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#parameter estimates
summary(pgls_eco1)
## 
## Call:
## pgls(formula = log10(eye_av) ~ log10(length_av) + Benthopelagic, 
##     data = species.comp, lambda = 1, param.CI = 0.95, bounds = list(lambda = c(1, 
##         1)))
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.092420 -0.049981  0.009349  0.043498  0.080494 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [Fix]  : 1.000
## delta  [Fix]  : 1.000
## 
## Coefficients:
##                   Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)      -1.287737   0.151661 -8.4909 1.158e-06 ***
## log10(length_av)  0.820905   0.098088  8.3691 1.360e-06 ***
## Benthopelagicyes  0.020845   0.036592  0.5697    0.5786    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.06001 on 13 degrees of freedom
## Multiple R-squared: 0.8448,  Adjusted R-squared: 0.821 
## F-statistic: 35.39 on 2 and 13 DF,  p-value: 5.496e-06

When lambda = 1, habitat has no effect on relative eye size.

Habitat (benthopelagic/pelagic) does not have a significant effect on relative eye size when lambda is set to 0 or 1. Thus, our finding is robust whatever the true value of lambda is. However, note that we only had n = 3 benthopelagic species in our dataset.